5 Integrating Minimalistic Localization and Navigation for People with

合集下载

晶体结构缺陷的类型

二按缺陷产生旳原因分类
晶体缺陷
辐照缺陷杂质缺陷
电荷缺陷热缺陷非化学计量缺陷
1. 热缺陷
定义:热缺陷亦称为本征缺陷，是指由热起伏旳原因所产生旳空位或间隙质点（原子或离子）。
类型:弗仑克尔缺陷（Frenkel defect）和肖特基缺陷（Schottky defect）
T E 热起伏(涨落) 原子脱离其平衡位置
面缺陷旳取向及分布与材料旳断裂韧性有关。
面缺陷－晶界
晶界示意图
亚晶界示意图
晶界: 晶界是两相邻晶粒间旳过渡界面。因为相邻晶粒间彼此位向各不相同，故晶界处旳原子排列与晶内不同，它们因同步受到相邻两侧晶粒不同位向旳综合影响，而做无规则排列或近似于两者取向旳折衷位置旳排列，这就形成了晶体中旳主要旳面缺陷。
-"extra" atoms positioned between atomic sites.
distortion of planes
selfinterstitiallids
Two outcomes if impurity (B) added to host (A):
• Solid solution of B in A (i.e., random dist. of point defects)
OR
Substitutional alloy (e.g., Cu in Ni)
Interstitial alloy (e.g., C in Fe)
Impurities in Ceramics
本章主要内容：
§2.1 晶体构造缺陷旳类型 §2. 2 点缺陷 §2.3 线缺陷 §2.4 面缺陷 §2.5 固溶体 §2.6 非化学计量化合物

Modeling the Spatial Dynamics of Regional Land Use_The CLUE-S Model

Modeling the Spatial Dynamics of Regional Land Use:The CLUE-S ModelPETER H.VERBURG*Department of Environmental Sciences Wageningen UniversityP.O.Box376700AA Wageningen,The NetherlandsandFaculty of Geographical SciencesUtrecht UniversityP.O.Box801153508TC Utrecht,The NetherlandsWELMOED SOEPBOERA.VELDKAMPDepartment of Environmental Sciences Wageningen UniversityP.O.Box376700AA Wageningen,The NetherlandsRAMIL LIMPIADAVICTORIA ESPALDONSchool of Environmental Science and Management University of the Philippines Los Ban˜osCollege,Laguna4031,Philippines SHARIFAH S.A.MASTURADepartment of GeographyUniversiti Kebangsaan Malaysia43600BangiSelangor,MalaysiaABSTRACT/Land-use change models are important tools for integrated environmental management.Through scenario analysis they can help to identify near-future critical locations in the face of environmental change.A dynamic,spatially ex-plicit,land-use change model is presented for the regional scale:CLUE-S.The model is speciﬁcally developed for the analysis of land use in small regions(e.g.,a watershed or province)at aﬁne spatial resolution.The model structure is based on systems theory to allow the integrated analysis of land-use change in relation to socio-economic and biophysi-cal driving factors.The model explicitly addresses the hierar-chical organization of land use systems,spatial connectivity between locations and stability.Stability is incorporated by a set of variables that deﬁne the relative elasticity of the actual land-use type to conversion.The user can specify these set-tings based on expert knowledge or survey data.Two appli-cations of the model in the Philippines and Malaysia are used to illustrate the functioning of the model and its validation.Land-use change is central to environmental man-agement through its inﬂuence on biodiversity,water and radiation budgets,trace gas emissions,carbon cy-cling,and livelihoods(Lambin and others2000a, Turner1994).Land-use planning attempts to inﬂuence the land-use change dynamics so that land-use conﬁg-urations are achieved that balance environmental and stakeholder needs.Environmental management and land-use planning therefore need information about the dynamics of land use.Models can help to understand these dynamics and project near future land-use trajectories in order to target management decisions(Schoonenboom1995).Environmental management,and land-use planning speciﬁcally,take place at different spatial and organisa-tional levels,often corresponding with either eco-re-gional or administrative units,such as the national or provincial level.The information needed and the man-agement decisions made are different for the different levels of analysis.At the national level it is often sufﬁ-cient to identify regions that qualify as“hot-spots”of land-use change,i.e.,areas that are likely to be faced with rapid land use conversions.Once these hot-spots are identiﬁed a more detailed land use change analysis is often needed at the regional level.At the regional level,the effects of land-use change on natural resources can be determined by a combina-tion of land use change analysis and speciﬁc models to assess the impact on natural resources.Examples of this type of model are water balance models(Schulze 2000),nutrient balance models(Priess and Koning 2001,Smaling and Fresco1993)and erosion/sedimen-tation models(Schoorl and Veldkamp2000).Most of-KEY WORDS:Land-use change;Modeling;Systems approach;Sce-nario analysis;Natural resources management*Author to whom correspondence should be addressed;email:pverburg@gissrv.iend.wau.nlDOI:10.1007/s00267-002-2630-x Environmental Management Vol.30,No.3,pp.391–405©2002Springer-Verlag New York Inc.ten these models need high-resolution data for land use to appropriately simulate the processes involved.Land-Use Change ModelsThe rising awareness of the need for spatially-ex-plicit land-use models within the Land-Use and Land-Cover Change research community(LUCC;Lambin and others2000a,Turner and others1995)has led to the development of a wide range of land-use change models.Whereas most models were originally devel-oped for deforestation(reviews by Kaimowitz and An-gelsen1998,Lambin1997)more recent efforts also address other land use conversions such as urbaniza-tion and agricultural intensiﬁcation(Brown and others 2000,Engelen and others1995,Hilferink and Rietveld 1999,Lambin and others2000b).Spatially explicit ap-proaches are often based on cellular automata that simulate land use change as a function of land use in the neighborhood and a set of user-speciﬁed relations with driving factors(Balzter and others1998,Candau 2000,Engelen and others1995,Wu1998).The speci-ﬁcation of the neighborhood functions and transition rules is done either based on the user’s expert knowl-edge,which can be a problematic process due to a lack of quantitative understanding,or on empirical rela-tions between land use and driving factors(e.g.,Pi-janowski and others2000,Pontius and others2000).A probability surface,based on either logistic regression or neural network analysis of historic conversions,is made for future conversions.Projections of change are based on applying a cut-off value to this probability sur-face.Although appropriate for short-term projections,if the trend in land-use change continues,this methodology is incapable of projecting changes when the demands for different land-use types change,leading to a discontinua-tion of the trends.Moreover,these models are usually capable of simulating the conversion of one land-use type only(e.g.deforestation)because they do not address competition between land-use types explicitly.The CLUE Modeling FrameworkThe Conversion of Land Use and its Effects(CLUE) modeling framework(Veldkamp and Fresco1996,Ver-burg and others1999a)was developed to simulate land-use change using empirically quantiﬁed relations be-tween land use and its driving factors in combination with dynamic modeling.In contrast to most empirical models,it is possible to simulate multiple land-use types simultaneously through the dynamic simulation of competition between land-use types.This model was developed for the national and con-tinental level,applications are available for Central America(Kok and Winograd2001),Ecuador(de Kon-ing and others1999),China(Verburg and others 2000),and Java,Indonesia(Verburg and others 1999b).For study areas with such a large extent the spatial resolution of analysis was coarse(pixel size vary-ing between7ϫ7and32ϫ32km).This is a conse-quence of the impossibility to acquire data for land use and all driving factors atﬁner spatial resolutions.A coarse spatial resolution requires a different data rep-resentation than the common representation for data with aﬁne spatial resolution.Inﬁne resolution grid-based approaches land use is deﬁned by the most dom-inant land-use type within the pixel.However,such a data representation would lead to large biases in the land-use distribution as some class proportions will di-minish and other will increase with scale depending on the spatial and probability distributions of the cover types(Moody and Woodcock1994).In the applications of the CLUE model at the national or continental level we have,therefore,represented land use by designating the relative cover of each land-use type in each pixel, e.g.a pixel can contain30%cultivated land,40%grass-land,and30%forest.This data representation is di-rectly related to the information contained in the cen-sus data that underlie the applications.For each administrative unit,census data denote the number of hectares devoted to different land-use types.When studying areas with a relatively small spatial ex-tent,we often base our land-use data on land-use maps or remote sensing images that denote land-use types respec-tively by homogeneous polygons or classiﬁed pixels. When converted to a raster format this results in only one, dominant,land-use type occupying one unit of analysis. The validity of this data representation depends on the patchiness of the landscape and the pixel size chosen. Most sub-national land use studies use this representation of land use with pixel sizes varying between a few meters up to about1ϫ1km.The two different data represen-tations are shown in Figure1.Because of the differences in data representation and other features that are typical for regional appli-cations,the CLUE model can not directly be applied at the regional scale.This paper describes the mod-iﬁed modeling approach for regional applications of the model,now called CLUE-S(the Conversion of Land Use and its Effects at Small regional extent). The next section describes the theories underlying the development of the model after which it is de-scribed how these concepts are incorporated in the simulation model.The functioning of the model is illustrated for two case-studies and is followed by a general discussion.392P.H.Verburg and othersCharacteristics of Land-Use SystemsThis section lists the main concepts and theories that are prevalent for describing the dynamics of land-use change being relevant for the development of land-use change models.Land-use systems are complex and operate at the interface of multiple social and ecological systems.The similarities between land use,social,and ecological systems allow us to use concepts that have proven to be useful for studying and simulating ecological systems in our analysis of land-use change (Loucks 1977,Adger 1999,Holling and Sanderson 1996).Among those con-cepts,connectivity is important.The concept of con-nectivity acknowledges that locations that are at a cer-tain distance are related to each other (Green 1994).Connectivity can be a direct result of biophysical pro-cesses,e.g.,sedimentation in the lowlands is a direct result of erosion in the uplands,but more often it is due to the movement of species or humans through the nd degradation at a certain location will trigger farmers to clear land at a new location.Thus,changes in land use at this new location are related to the land-use conditions in the other location.In other instances more complex relations exist that are rooted in the social and economic organization of the system.The hierarchical structure of social organization causes some lower level processes to be constrained by higher level dynamics,e.g.,the establishments of a new fruit-tree plantation in an area near to the market might in ﬂuence prices in such a way that it is no longer pro ﬁtable for farmers to produce fruits in more distant areas.For studying this situation an-other concept from ecology,hierarchy theory,is use-ful (Allen and Starr 1982,O ’Neill and others 1986).This theory states that higher level processes con-strain lower level processes whereas the higher level processes might emerge from lower level dynamics.This makes the analysis of the land-use system at different levels of analysis necessary.Connectivity implies that we cannot understand land use at a certain location by solely studying the site characteristics of that location.The situation atneigh-Figure 1.Data representation and land-use model used for respectively case-studies with a national/continental extent and local/regional extent.Modeling Regional Land-Use Change393boring or even more distant locations can be as impor-tant as the conditions at the location itself.Land-use and land-cover change are the result of many interacting processes.Each of these processes operates over a range of scales in space and time.These processes are driven by one or more of these variables that inﬂuence the actions of the agents of land-use and cover change involved.These variables are often re-ferred to as underlying driving forces which underpin the proximate causes of land-use change,such as wood extraction or agricultural expansion(Geist and Lambin 2001).These driving factors include demographic fac-tors(e.g.,population pressure),economic factors(e.g., economic growth),technological factors,policy and institutional factors,cultural factors,and biophysical factors(Turner and others1995,Kaimowitz and An-gelsen1998).These factors inﬂuence land-use change in different ways.Some of these factors directly inﬂu-ence the rate and quantity of land-use change,e.g.the amount of forest cleared by new incoming migrants. Other factors determine the location of land-use change,e.g.the suitability of the soils for agricultural land use.Especially the biophysical factors do pose constraints to land-use change at certain locations, leading to spatially differentiated pathways of change.It is not possible to classify all factors in groups that either inﬂuence the rate or location of land-use change.In some cases the same driving factor has both an inﬂu-ence on the quantity of land-use change as well as on the location of land-use change.Population pressure is often an important driving factor of land-use conver-sions(Rudel and Roper1997).At the same time it is the relative population pressure that determines which land-use changes are taking place at a certain location. Intensively cultivated arable lands are commonly situ-ated at a limited distance from the villages while more extensively managed grasslands are often found at a larger distance from population concentrations,a rela-tion that can be explained by labor intensity,transport costs,and the quality of the products(Von Thu¨nen 1966).The determination of the driving factors of land use changes is often problematic and an issue of dis-cussion(Lambin and others2001).There is no unify-ing theory that includes all processes relevant to land-use change.Reviews of case studies show that it is not possible to simply relate land-use change to population growth,poverty,and infrastructure.Rather,the inter-play of several proximate as well as underlying factors drive land-use change in a synergetic way with large variations caused by location speciﬁc conditions (Lambin and others2001,Geist and Lambin2001).In regional modeling we often need to rely on poor data describing this complexity.Instead of using the under-lying driving factors it is needed to use proximate vari-ables that can represent the underlying driving factors. Especially for factors that are important in determining the location of change it is essential that the factor can be mapped quantitatively,representing its spatial vari-ation.The causality between the underlying driving factors and the(proximate)factors used in modeling (in this paper,also referred to as“driving factors”) should be certiﬁed.Other system properties that are relevant for land-use systems are stability and resilience,concepts often used to describe ecological systems and,to some extent, social systems(Adger2000,Holling1973,Levin and others1998).Resilience refers to the buffer capacity or the ability of the ecosystem or society to absorb pertur-bations,or the magnitude of disturbance that can be absorbed before a system changes its structure by changing the variables and processes that control be-havior(Holling1992).Stability and resilience are con-cepts that can also be used to describe the dynamics of land-use systems,that inherit these characteristics from both ecological and social systems.Due to stability and resilience of the system disturbances and external in-ﬂuences will,mostly,not directly change the landscape structure(Conway1985).After a natural disaster lands might be abandoned and the population might tempo-rally migrate.However,people will in most cases return after some time and continue land-use management practices as before,recovering the land-use structure (Kok and others2002).Stability in the land-use struc-ture is also a result of the social,economic,and insti-tutional structure.Instead of a direct change in the land-use structure upon a fall in prices of a certain product,farmers will wait a few years,depending on the investments made,before they change their cropping system.These characteristics of land-use systems provide a number requirements for the modelling of land-use change that have been used in the development of the CLUE-S model,including:●Models should not analyze land use at a single scale,but rather include multiple,interconnected spatial scales because of the hierarchical organization of land-use systems.●Special attention should be given to the drivingfactors of land-use change,distinguishing drivers that determine the quantity of change from drivers of the location of change.●Sudden changes in driving factors should not di-rectly change the structure of the land-use system asa consequence of the resilience and stability of theland-use system.394P.H.Verburg and others●The model structure should allow spatial interac-tions between locations and feedbacks from higher levels of organization.Model DescriptionModel StructureThe model is sub-divided into two distinct modules,namely a non-spatial demand module and a spatially explicit allocation procedure (Figure 2).The non-spa-tial module calculates the area change for all land-use types at the aggregate level.Within the second part of the model these demands are translated into land-use changes at different locations within the study region using a raster-based system.For the land-use demand module,different alterna-tive model speci ﬁcations are possible,ranging from simple trend extrapolations to complex economic mod-els.The choice for a speci ﬁc model is very much de-pendent on the nature of the most important land-use conversions taking place within the study area and the scenarios that need to be considered.Therefore,the demand calculations will differ between applications and scenarios and need to be decided by the user for the speci ﬁc situation.The results from the demandmodule need to specify,on a yearly basis,the area covered by the different land-use types,which is a direct input for the allocation module.The rest of this paper focuses on the procedure to allocate these demands to land-use conversions at speci ﬁc locations within the study area.The allocation is based upon a combination of em-pirical,spatial analysis,and dynamic modelling.Figure 3gives an overview of the procedure.The empirical analysis unravels the relations between the spatial dis-tribution of land use and a series of factors that are drivers and constraints of land use.The results of this empirical analysis are used within the model when sim-ulating the competition between land-use types for a speci ﬁc location.In addition,a set of decision rules is speci ﬁed by the user to restrict the conversions that can take place based on the actual land-use pattern.The different components of the procedure are now dis-cussed in more detail.Spatial AnalysisThe pattern of land use,as it can be observed from an airplane window or through remotely sensed im-ages,reveals the spatial organization of land use in relation to the underlying biophysical andsocio-eco-Figure 2.Overview of the modelingprocedure.Figure 3.Schematic represen-tation of the procedure to allo-cate changes in land use to a raster based map.Modeling Regional Land-Use Change395nomic conditions.These observations can be formal-ized by overlaying this land-use pattern with maps de-picting the variability in biophysical and socio-economic conditions.Geographical Information Systems(GIS)are used to process all spatial data and convert these into a regular grid.Apart from land use, data are gathered that represent the assumed driving forces of land use in the study area.The list of assumed driving forces is based on prevalent theories on driving factors of land-use change(Lambin and others2001, Kaimowitz and Angelsen1998,Turner and others 1993)and knowledge of the conditions in the study area.Data can originate from remote sensing(e.g., land use),secondary statistics(e.g.,population distri-bution),maps(e.g.,soil),and other sources.To allow a straightforward analysis,the data are converted into a grid based system with a cell size that depends on the resolution of the available data.This often involves the aggregation of one or more layers of thematic data,e.g. it does not make sense to use a30-m resolution if that is available for land-use data only,while the digital elevation model has a resolution of500m.Therefore, all data are aggregated to the same resolution that best represents the quality and resolution of the data.The relations between land use and its driving fac-tors are thereafter evaluated using stepwise logistic re-gression.Logistic regression is an often used method-ology in land-use change research(Geoghegan and others2001,Serneels and Lambin2001).In this study we use logistic regression to indicate the probability of a certain grid cell to be devoted to a land-use type given a set of driving factors following:LogͩP i1ϪP i ͪϭ␤0ϩ␤1X1,iϩ␤2X2,i......ϩ␤n X n,iwhere P i is the probability of a grid cell for the occur-rence of the considered land-use type and the X’s are the driving factors.The stepwise procedure is used to help us select the relevant driving factors from a larger set of factors that are assumed to inﬂuence the land-use pattern.Variables that have no signiﬁcant contribution to the explanation of the land-use pattern are excluded from theﬁnal regression equation.Where in ordinal least squares regression the R2 gives a measure of modelﬁt,there is no equivalent for logistic regression.Instead,the goodness ofﬁt can be evaluated with the ROC method(Pontius and Schnei-der2000,Swets1986)which evaluates the predicted probabilities by comparing them with the observed val-ues over the whole domain of predicted probabilities instead of only evaluating the percentage of correctly classiﬁed observations at aﬁxed cut-off value.This is an appropriate methodology for our application,because we will use a wide range of probabilities within the model calculations.The inﬂuence of spatial autocorrelation on the re-gression results can be minimized by only performing the regression on a random sample of pixels at a certain minimum distance from one another.Such a selection method is adopted in order to maximize the distance between the selected pixels to attenuate the problem associated with spatial autocorrelation.For case-studies where autocorrelation has an important inﬂuence on the land-use structure it is possible to further exploit it by incorporating an autoregressive term in the regres-sion equation(Overmars and others2002).Based upon the regression results a probability map can be calculated for each land-use type.A new probabil-ity map is calculated every year with updated values for the driving factors that are projected to change in time,such as the population distribution or accessibility.Decision RulesLand-use type or location speciﬁc decision rules can be speciﬁed by the user.Location speciﬁc decision rules include the delineation of protected areas such as nature reserves.If a protected area is speciﬁed,no changes are allowed within this area.For each land-use type decision rules determine the conditions under which the land-use type is allowed to change in the next time step.These decision rules are implemented to give certain land-use types a certain resistance to change in order to generate the stability in the land-use structure that is typical for many landscapes.Three different situations can be distinguished and for each land-use type the user should specify which situation is most relevant for that land-use type:1.For some land-use types it is very unlikely that theyare converted into another land-use type after their ﬁrst conversion;as soon as an agricultural area is urbanized it is not expected to return to agriculture or to be converted into forest cover.Unless a de-crease in area demand for this land-use type occurs the locations covered by this land use are no longer evaluated for potential land-use changes.If this situation is selected it also holds that if the demand for this land-use type decreases,there is no possi-bility for expansion in other areas.In other words, when this setting is applied to forest cover and deforestation needs to be allocated,it is impossible to reforest other areas at the same time.2.Other land-use types are converted more easily.Aswidden agriculture system is most likely to be con-verted into another land-use type soon after its396P.H.Verburg and othersinitial conversion.When this situation is selected for a land-use type no restrictions to change are considered in the allocation module.3.There is also a number of land-use types that oper-ate in between these two extremes.Permanent ag-riculture and plantations require an investment for their establishment.It is therefore not very likely that they will be converted very soon after into another land-use type.However,in the end,when another land-use type becomes more pro ﬁtable,a conversion is possible.This situation is dealt with by de ﬁning the relative elasticity for change (ELAS u )for the land-use type into any other land use type.The relative elasticity ranges between 0(similar to Situation 2)and 1(similar to Situation 1).The higher the de ﬁned elasticity,the more dif ﬁcult it gets to convert this land-use type.The elasticity should be de ﬁned based on the user ’s knowledge of the situation,but can also be tuned during the calibration of the petition and Actual Allocation of Change Allocation of land-use change is made in an iterative procedure given the probability maps,the decision rules in combination with the actual land-use map,and the demand for the different land-use types (Figure 4).The following steps are followed in the calculation:1.The ﬁrst step includes the determination of all grid cells that are allowed to change.Grid cells that are either part of a protected area or under a land-use type that is not allowed to change (Situation 1,above)are excluded from further calculation.2.For each grid cell i the total probability (TPROP i,u )is calculated for each of the land-use types u accord-ing to:TPROP i,u ϭP i,u ϩELAS u ϩITER u ,where ITER u is an iteration variable that is speci ﬁc to the land use.ELAS u is the relative elasticity for change speci ﬁed in the decision rules (Situation 3de-scribed above)and is only given a value if grid-cell i is already under land use type u in the year considered.ELAS u equals zero if all changes are allowed (Situation 2).3.A preliminary allocation is made with an equalvalue of the iteration variable (ITER u )for all land-use types by allocating the land-use type with the highest total probability for the considered grid cell.This will cause a number of grid cells to change land use.4.The total allocated area of each land use is nowcompared to the demand.For land-use types where the allocated area is smaller than the demanded area the value of the iteration variable is increased.For land-use types for which too much is allocated the value is decreased.5.Steps 2to 4are repeated as long as the demandsare not correctly allocated.When allocation equals demand the ﬁnal map is saved and the calculations can continue for the next yearly timestep.Figure 5shows the development of the iteration parameter ITER u for different land-use types during asimulation.Figure 4.Representation of the iterative procedure for land-use changeallocation.Figure 5.Change in the iteration parameter (ITER u )during the simulation within one time-step.The different lines rep-resent the iteration parameter for different land-use types.The parameter is changed for all land-use types synchronously until the allocated land use equals the demand.Modeling Regional Land-Use Change397Multi-Scale CharacteristicsOne of the requirements for land-use change mod-els are multi-scale characteristics.The above described model structure incorporates different types of scale interactions.Within the iterative procedure there is a continuous interaction between macro-scale demands and local land-use suitability as determined by the re-gression equations.When the demand changes,the iterative procedure will cause the land-use types for which demand increased to have a higher competitive capacity (higher value for ITER u )to ensure enough allocation of this land-use type.Instead of only being determined by the local conditions,captured by the logistic regressions,it is also the regional demand that affects the actually allocated changes.This allows the model to “overrule ”the local suitability,it is not always the land-use type with the highest probability according to the logistic regression equation (P i,u )that the grid cell is allocated to.Apart from these two distinct levels of analysis there are also driving forces that operate over a certain dis-tance instead of being locally important.Applying a neighborhood function that is able to represent the regional in ﬂuence of the data incorporates this type of variable.Population pressure is an example of such a variable:often the in ﬂuence of population acts over a certain distance.Therefore,it is not the exact location of peoples houses that determines the land-use pattern.The average population density over a larger area is often a more appropriate variable.Such a population density surface can be created by a neighborhood func-tion using detailed spatial data.The data generated this way can be included in the spatial analysis as anotherindependent factor.In the application of the model in the Philippines,described hereafter,we applied a 5ϫ5focal ﬁlter to the population map to generate a map representing the general population pressure.Instead of using these variables,generated by neighborhood analysis,it is also possible to use the more advanced technique of multi-level statistics (Goldstein 1995),which enable a model to include higher-level variables in a straightforward manner within the regression equa-tion (Polsky and Easterling 2001).Application of the ModelIn this paper,two examples of applications of the model are provided to illustrate its function.TheseTable nd-use classes and driving factors evaluated for Sibuyan IslandLand-use classes Driving factors (location)Forest Altitude (m)GrasslandSlope Coconut plantation AspectRice ﬁeldsDistance to town Others (incl.mangrove and settlements)Distance to stream Distance to road Distance to coast Distance to port Erosion vulnerability GeologyPopulation density(neighborhood 5ϫ5)Figure 6.Location of the case-study areas.398P.H.Verburg and others。

Integrating Time Alignment and Self-Organizing Maps for Classifying Curves

Integrating Time Alignment and Self-Organizing Mapsfor Classifying CurvesElvira Romano and Germana ScepiDipartimento di Matematica e Statistica – Università “Federico II” di NapoliVia Cintia, Monte S. Angelo – 80126 NapoliKeywords: Classification, functional data, time series, dissimilarity.1.IntroductionClustering time series has become in recent years a topic of great interest in a wide range of fields. The several approaches differ mainly in their notion of similarity (for a review see Focardi, 2001). Most researches use the Euclidean distance or some variation of it because of its easy implementation, even if it is very sensitive to temporal axis alignment.Furthermore, there are many applications where it is demonstrated that the Euclidean distances between raw data fail to capture the notion of similarity. The principal reason why Euclidean distance may fail to produce an intuitively correct measure of similarity between two sequences is that it is very sensitive to small distortions in the time axis as, for example, in the case of two sequences having approximately the same overall shape but not aligned in time axis.A method that allows this elastic shifting of the X-axis is desired in order to detect similar shapes with different phases. For this purpose, the Dynamical Time Warping (DTW) distance has been recently introduced (Berndt, Clifford, 1994), technique that was already known in the speech processing community (Sakoe, Chiba, 1978; Rabiner, Juang, 1993).Nevertheless the DTW algorithm can produce incorrect results in presence of salient features or noise in the data and the algorithm’s time complexity causes a problem in a way that “…performance on very large databases may be a limitation”.Morlini et al. (2005) proposes a modification of this algorithm that considers a smoothed version of the data and demonstrate that their approach allows to obtain points which are less noisy and dependent on the overall shape of the series. The clustering algorithms proposed in this approach are hierarchical clustering and K-means algorithms.The current paper proposes a new approach based on the implementation of the DTW distance in a Self Organizing Map algorithm (Kohonen, 2001) with the aim of classifying a set of curves. To show the results of this approach, we illustrate an application of our method on simulated data; while in the extended paper version we will propose an application on topographic real data.2.A new approach for classifying curvesSuppose we have several time series. Let us consider, for example Q and C, two time series of length n and m respectively:Q = q 1,q 2,…,q i ,…,q nC = c 1,c 2,…,c j ,…,c mThe first step of our approach consists in smoothing each series by a piecewise linear or cubic spline. Therefore our starting data are a set of curves, in the example:Q ’ = ''''12,,...,,...,i n q q q q C ’ =''''12 ,,...,,...,j m c c c cTo align the two obtained sequences using DTW, we construct an n-by-m matrix where the (i th ,j th ) element of the matrix contains the Euclidean distance between thetwo points and .Each matrix element (i,j ) corresponds to the alignment between thepoints.''i j (,)d q c 'i q 'j c A warping path, W , is a contiguous set of matrix elements that defines a mapping between Q ’ and C ’. The k -th element of W is defined as w k = (i,j )k , so we have:12,,...,,..., max(,)-1k K W w w w w m n K m n =≤≤+The warping path is typically subjected to several constraints-Boundary conditions: w 1 = (1,1) and w K = (m ,n ). Simply stated, this requires the warping path to start and finish in diagonally opposite corner cells of the matrix.- Continuity: Given w k = (a ,b ) then w k-1 = (a’,b’), where a–a' ≤1 and b-b' ≤1. This restricts the allowable steps in the warping path to adjacent cells (including diagonally adjacent cells).- Monotonicity: Given w k = (a ,b ) then w k-1 = (a',b'), where a–a' ≥0 and b-b'≥ 0. This forces the points in W to be monotonically spaced in time.We are interested only in the path that minimizes the warping cost:DTWC(',')min Q C = (1)Therefore in our approach, the data are the smoothed values of sequences and the dissimilarity between two elements is the Dynamic Time Warping Cost (DTWC). The clustering method is based on an adaptation of the Kohonen’s SOM algorithm for dissimilarity data (Golli et al, 2004).The SOM algorithm consists of neurons organized on a regular low dimensional map. More formally, the map is described by a graph (N ,Γ). N is a set of interconnectedneurons having a discrete topology defined by Γ. For each pair of neurons on the map, the distance is defined as the shortest path between them on the graph. This distance imposes a neighbourhood relation between neurons.The Dissimilarity SOM algorithm (DSOM) is an adaptation of the Kohonen’s SOM algorithm for dissimilarity data. It is a batch iterative algorithm in which the whole data set (Ω) is initially presented on the map. We denote with z l (l=1,…,N) the generic element of Ω and z l is the representation of this element in representative space D on which dissimilarity (denoted d ) is defined.Each neuron x is represented by a set of M elements of Ω , m 1,…,m g ,…,m M , called prototypes, where m g is a vector of z l element.In DSOM the prototypes associated to neurons as well as the neighbourhood function are evolving with the iterations. It starts by an initialization phase, in which the value of M is randomly chosen.The algorithm alternates affectation phases and representation phases until convergence. In the first phase each initial observation is assigned to the winning prototype according to the following assignment function:()(arg min ,T l g M)l g f z d z ∈=m (2)where the adequacy function is:()()()2,,s r T T l g l s g r r M z m d z m K d ϑ∈∈=∑∑,z z (3)with (),T g r K ϑ is the neighbourhood kernel around the neuron r and ,l s z z are the representations of the elements in the space D .At the generic h -th iteration we assign an observation to the winning prototype with the (2) and define the cluster associated to this prototype at the iteration h . The main drawback of the DSOM algorithm is the cost induced by the representation phase. A fast version of the DSOM algorithm that allows a an important reduction of its theoretical cost has been proposed by Conan-Guez et al. (2005).In our approach we aim to classify a set of curves by using the described clustering algorithm and the DTWC. Therefore in our algorithm the smoothed time series are classified by substituting the distance d in (3) with the DTWC (1).This approach allows us to have an easy visualization of data and it is computationally more efficient of the classical clustering algorithms, it deals with time series drawn from large data sets. The visualization of time series is very important for the detection of their own characteristics and gives us some information for representing each class.3. Experimental ResultsFor a first evaluation of our approach, we propose a simulation study on a small data set of 130 time series. We have generated 130 time series (Fig.1) of length 100 and in particular i) 60 time series with increasing trend, ii) 30 time series with a seasonal component only and iii) 40 time series with decreasing trend .Fig.1 The simulated time seriesWithout warping, the k-means algorithm is not often able to distinguish class i) from class ii), with a general misclassification of 53%.We wrote a Matlab program for generating the time series, smoothing each series with a cubic spline and implementing the DTW algorithm. Finally, we clustered the series on the basis of the DTWC with the DSOM algorithm.We have repeated the simulation study 150 times with different values of the smoothing parameter λ ranging from 0.05 to 0.20.The results show that (Fig.2) only the 10% of time series are misclassified (for λ=0.1) and there are very few cases of confusion between class i) and class ii).Fig.2 The classification results with the smoothed time series4. ConclusionsIn this short version of the paper, we have introduced a new approach for clustering smoothed time series, based on the joint use of the Dynamic Time Warping distance and the Dissimilarity SOM algorithm.This algorithm seems particularly promising in data mining problems and it can be applied on not aligned time points with a good visualisation of results.The forthcoming paper includes a more detailed analysis of our approach and, in particular, an application on a large set of real data, which is needed to investigate the robustness of the proposed approach in presence of irregular sampled data. A comparison, on the same data, with the algorithm proposed by Morlini et al. (2005) will be performed.In the further researches we aim to define a non parametric model for characterizing each obtained cluster. In other words, we will search for each cluster of smoothed time series a non parametric function synthesizing its elements.Main ReferencesBerndt, D., Clifford, J. (1994). Using dynamic time warping to find patterns in time series AAI -94 Workshop on Knowledge Discovery in Databases,229–248.Conan-Guez B., Rossi F., Golli A.E. (2005). A Fast Algorithm for the Self-Organizing Map on Dissimilarity Data,WSOM’05 Proceedings.Focardi S.M. (2001). Clustering delle serie storiche economiche: applicazioni e questioni computazionali, Technical Report, Supercalcolo in Economia e in Finanza Milano. Golli A.E., Conan-Guez B, Rossi F. (2004). Self-organizing maps and symbolic data CLUEB, Journal of Symbolic Data Analysis, 2, n.1, ISSN 1723–5081.Kohonen T. (2001). Self-Organizing Maps, Springer Series in Information Sciences, Springer.Molini I. (2005). On the Dynamic Time Warping for Computing the Dissimilarity Between Curves, Vichi et al. eds., New Developments in Classification and Data Analysis, Proceedings of the Meeting of the Classification and Data Analysis Group, Università di Bologna, Settembre 2003.Rabiner, L., Juang, B. (1993). Fundamentals of speech recognition, Englewood Cliffs, N.J, Prentice Hall.Sakoe, H., Chiba, S. (1978). Dynamic programming algorithm optimization for spoken word recognition. IEEE Trans. Acoustics, Speech, and Signal Processing., Volume 26, pp 143-165.。

PREFACE

2
Nara Women’s University
3
Table of Contents
Session 1 Session 2 Session 3 Session 4 Session 5 Session 6 Session 7 Session 8 Session 9 Session 10 Session 11 Session 12 Session 13 Session 14
General Session Young Researchers Session Approximate Algebraic Computation Computational Algebraic Structures and Engineering Applications Computer Algebra and Coding Theory Computer Algebra in Quantum Information and Computation Computer Algebra in the Biological Sciences Computational Topology and Geometry Computer Algebra in Education Handling Large Expressions in Symbolic Computation High-Performance Computer Algebra Newton and Hensel Techniques in Scientiﬁc Computing Parametric and Nonconvex Constraint Solving Pen-Based Mathematical Computing
Abstracts of Presentations

集成梯度特征归属方法-概述说明以及解释

集成梯度特征归属方法-概述说明以及解释1.引言1.1 概述在概述部分，你可以从以下角度来描述集成梯度特征归属方法的背景和重要性：集成梯度特征归属方法是一种用于分析和解释机器学习模型预测结果的技术。

随着机器学习的快速发展和广泛应用，对于模型的解释性需求也越来越高。

传统的机器学习模型通常被认为是“黑盒子”，即无法解释模型做出预测的原因。

这限制了模型在一些关键应用领域的应用，如金融风险评估、医疗诊断和自动驾驶等。

为了解决这个问题，研究人员提出了各种机器学习模型的解释方法，其中集成梯度特征归属方法是一种非常受关注和有效的技术。

集成梯度特征归属方法能够为机器学习模型的预测结果提供可解释的解释，从而揭示模型对于不同特征的关注程度和影响力。

通过分析模型中每个特征的梯度值，可以确定该特征在预测中扮演的角色和贡献度，从而帮助用户理解模型的决策过程。

这对于模型的评估、优化和改进具有重要意义。

集成梯度特征归属方法的应用广泛，不仅适用于传统的机器学习模型，如决策树、支持向量机和逻辑回归等，也可以应用于深度学习模型，如神经网络和卷积神经网络等。

它能够为各种类型的特征，包括数值型特征和类别型特征，提供有益的信息和解释。

本文将对集成梯度特征归属方法的原理、应用优势和未来发展进行详细阐述，旨在为读者提供全面的了解和使用指南。

在接下来的章节中，我们将首先介绍集成梯度特征归属方法的基本原理和算法，然后探讨应用该方法的优势和实际应用场景。

最后，我们将总结该方法的重要性，并展望未来该方法的发展前景。

1.2文章结构文章结构内容应包括以下内容：文章的结构部分主要是对整篇文章的框架进行概述，指导读者在阅读过程中能够清晰地了解文章的组织结构和内容安排。

第一部分是引言，介绍了整篇文章的背景和意义。

其中，1.1小节概述文章所要讨论的主题，简要介绍了集成梯度特征归属方法的基本概念和应用领域。

1.2小节重点在于介绍文章的结构，将列出本文各个部分的标题和内容概要，方便读者快速了解文章的大致内容。

基于声学指数的神农架国家公园声音多样性动态变化

果显示ＡＣＩ指数不能很好地反映日变化趋势，但ＢＩ指数和ＮＤＳＩ指数具有明显的日变化趋势，且变化趋势符合
物种黎明／黄昏合唱的习性；声学指数随海拔梯度的空间变化结果表明，ＡＣＩ、ＢＩ指数在中海拔区域具有最大值，
且ＡＣＩ指数与海拔相关性较强，ＮＤＳＩ指数没有显著的变化趋势。【结论】ＢＩ、ＮＤＳＩ指数能较好地反映动物声音
：【】ＡｂｓｔｒａｃｔＯｂｊｅｃｔｉｖｅＴｈｅｓｔｕｄｙａｉｍｓｔｏｅｖａｌｕａｔｅｔｈｅｒｅｓｐｏｎｓｅｏｆａｃｏｕｓｔｉｃｉｎｄｉｃｅｓｔｏｔｈｅｄｙｎａｍｉｃｃｈａｎｇｅｓｏｆａｎｉｍａｌ，ｓｏｕｎｄｄｉｖｅｒｓｉｔｙｆｕｒｔｈｅｒｔｏｅｘｐｌｏｒｅｔｈｅｃｈａｒａｃｔｅｒｉｓｔｉｃｓｏｆｔｈｅｖａｒｉａｔｉｏｎｏｆａｎｉｍａｌｓｏｕｎｄｄｉｖｅｒｓｉｔｙｉｎＳｈｅｎｎｏｎｇｊｉａＮａｔｉｏｎａｌ，，【】ＰａｒｋＣｈｉｎａｉｎｏｒｄｅｒｔｏｐｒｏｖｉｄｅａｑｕａｎｔｉｔａｔｉｖｅｂａｓｉｓｆｏｒｔｈｅｌｏｃａｌｅｃｏｌｏｇｉｃａｌｐｒｏｔｅｃｔｉｏｎ．ＭｅｔｈｏｄＷｅｄｅｐｌｏｙｅｄｎｉｎｅ，ｓｏｕｎｄｒｅｃｏｒｄｉｎｇｅｑｕｉｐｍｅｎｔｓｉｎｎｉｎｅｓａｍｐｌｉｎｇｓｉｔｅｓｉｎＳｈｅｎｎｏｎｇｊｉａＮａｔｉｏｎａｌＰａｒｋａｎｄｓｏｕｎｄｒｅｃｏｒｄｉｎｇｄａｔａｆｒｏｍＭａｙｔｏ（），Ｊｕｌｙ２０２１ｗｅｒｅｏｂｔａｉｎｅｄ．ＡｔｉｍｅｓｅｒｉｅｓｏｆｅｃｏａｃｏｓｕｔｉｃｉｎｄｉｃｅｓｉｎｃｌｕｄｉｎｇａｃｏｕｓｔｉｃｃｏｍｐｌｅｘｉｔｙｉｎｄｅｘＡＣＩｂｉｏａｃｏｕｓｔｉｃ（），（）ｉｎｄｅｘＢＩｎｏｒｍａｌｉｚｅｄｄｉｆｆｅｒｅｎｃｅｓｏｕｎｄｓｃａｐｅｉｎｄｅｘＮＤＳＩｗｅｒｅｅｘｔｒａｃｔｅｄｆｒｏｍｔｈｅｒｅｃｏｒｄｉｎｇｄａｔａａｆｔｅｒｎｏｉｓｅ

电容去离子除氯电极的构建及其脱盐性能研究进展

物理化学学报Acta Phys. -Chim. Sin. 2022, 38 (5), 2006037 (1 of 12)Received: June 12, 2020; Revised: July 2, 2020; Accepted: July 7, 2020; Published online: July 13, 2020. *Correspondingauthor.Email:**************.cn.The project was supported by the National Natural Science Foundation of China (21777118). 国家自然科学基金(21777118)资助项目© Editorial office of Acta Physico-Chimica Sinica[Review] doi: 10.3866/PKU.WHXB202006037 Research Progress in Chlorine Ion Removal Electrodes for Desalination by Capacitive DeionizationYuecheng Xiong 1, Fei Yu 2, Jie Ma 1,3,*1 Key Laboratory of Yangtze River Water Environment, Tongji University, Shanghai 200092, China.2 College of Marine Ecology and Environment, Shanghai Ocean University, Shanghai 201306, China.3 Shanghai Institute of Pollution Control and Ecological Security, Shanghai 200092, China .Abstract: Sustainable freshwater supply is a grave challenge to the society because of the severe water scarcity and global pollution. Seawater is an inexhaustible source of industrial and potable water. The relevant desalination technologies with a high market share include reverse osmosis and thermal distillation, which are energy-intensive. Capacitive deionization (CDI) is a desalination technology that is gaining extensive attention because of its low energy consumption and low chemical intensity. In CDI, charged species are removed from the aqueous environment via applying a voltage onto the anode and cathode. For desalination, Na + and Cl − ions are removed by the cathode and anode, respectively. With the boom in electrode materials for rechargeable batteries, the Na + removal electrode (cathode) hasevolved from a carbon-based electrode to a faradaic electrode, and the desalination performance of CDI has also been significantly enhanced. A conventional carbon-based electrode captures ions in the electrical double layer (EDL) and suffers from low charge efficiency, thus being unsuitable for use in water with high salinity. On the other hand, a faradaic electrode stores Na + ions through a reversible redox process or intercalation, leading to high desalination capacity.However, the Cl − removal electrode (anode) has not yet seen notable development. Most research groups employ activated carbon to remove Cl −, and therefore, summarizing Cl − storage electrodes for CDI is necessary to guide the design of electrode systems with better desalination performance. First, this review outlines the evolution of CDI configuration based on the electrode materials, suggesting that the anode and cathode are of equal importance in CDI. Second, a systematic summary of the anode materials used in CDI and a comparison of the characteristics of different electrodes, including those based on Ag/AgCl, Bi/BiOCl, 2-dimensional (2D) materials (layered double hydroxide (LDH) and MXene), redox polymers, and electrolytes, are presented. Then, the underlying mechanism for Cl − storage is refined. Similar to the case of Na + storage, traditional carbon electrodes store Cl- via electrosorption based on the EDL. Ag/AgCl and Bi/BiOCl remove Cl − through a conversion reaction, i.e ., phase transformation during the reaction with Cl −. 2D materials store Cl − in the space between adjacent layers, a process referred as ion intercalation, with layered double hydroxide (LDH) and MXene showing higher Cl − storage potential. Redox polymers and electrolytes allow for Cl − storage via redox reactions. Among all the materials mentioned above, Bi/BiOCl and LDH are the most promising for the construction of CDI anodes because of their high capacity and low cost. Finally, to spur the development of novel anodes for CDI, the electrodes applied in a chlorine ion battery are introduced. This is the first paper to comb through reports on the development of anode materials for CDI, thus laying the theoretical foundation for future materials design. Key Words: Capacitive deionization; Desalination; Anode; Chlorineion; Battery电容去离子除氯电极的构建及其脱盐性能研究进展熊岳城1，于飞2，马杰1,3,*1同济大学长江水环境教育部重点实验室，上海 2000922上海海洋大学海洋生态与环境学院，上海 2013063上海污染控制与生态安全研究院，上海 200092摘要：电容去离子技术(Capacitive deionization，CDI)是一种新兴的脱盐技术，通过在电极两端施加较低的外加电场除去水中的带电离子和分子，由于其较低的能耗和可持续性而备受关注。

Accurate Passive Location Estimation Using TOA Measurements

Accurate Passive Location Estimation Using TOA MeasurementsJunyang Shen,Andreas F.Molisch,Fellow,IEEE,and Jussi Salmi,Member,IEEEAbstract—Localization of objects is fast becoming a major aspect of wireless technologies,with applications in logistics, surveillance,and emergency response.Time-of-arrival(TOA) localization is ideally suited for high-precision localization of objects in particular in indoor environments,where GPS is not available.This paper considers the case where one transmitter and multiple,distributed,receivers are used to estimate the location of a passive(reﬂecting)object.It furthermore focuses on the situation when the transmitter and receivers can be synchronized,so that TOA(as opposed to time-difference-of-arrival(TDOA))information can be used.We propose a novel, Two-Step estimation(TSE)algorithm for the localization of the object.We then derive the Cramer-Rao Lower Bound(CRLB) for TOA and show that it is an order of magnitude lower than the CRLB of TDOA in typical setups.The TSE algorithm achieves the CRLB when the TOA measurements are subject to small Gaussian-distributed errors,which is veriﬁed by analytical and simulation results.Moreover,practical measurement results show that the estimation error variance of TSE can be33dB lower than that of TDOA based algorithms.Index Terms—TOA,TDOA,location estimation,CRLB.I.I NTRODUCTIONO BJECT location estimation has recently received inten-sive interests for a large variety of applications.For example,localization of people in smoke-ﬁlled buildings can be life-saving[1];positioning techniques also provide useful location information for search-and-rescue[2],logistics[3], and security applications such as localization of intruders[4].A variety of localization techniques have been proposed in the literature,which differ by the type of information and system parameters that are used.The three most important kinds utilize the received signal strength(RSS)[5],angle of arrival(AOA)[6],and signal propagation time[7],[8],[9], respectively.RSS algorithms use the received signal power for object positioning;their accuracies are limited by the fading of wireless signals[5].AOA algorithms require either directional antennas or receiver antenna arrays1.Signal-propagation-time based algorithms estimate the object location using the time it takes the signal to travel from the transmitter to the target and from there to the receivers.They achieve very accurate Manuscript received April15,2011;revised September28,2011and Jan-uary18,2012;accepted February12,2012.The associate editor coordinating the review of this paper and approving it for publication was X.Wang.J.Shen and A.F.Molisch are,and J.Salmi was with the Department of Electrical Engineering,Viterbi School of Engineering,University of Southern California(e-mail:{junyangs,molisch,salmi}@).J.Salmi is currently with Aalto University,SMARAD CoE,Espoo,Finland.This paper is partially supported by the Ofﬁce of Naval Research(ONR) under grant10599363.Part of this work was presented in the IEEE Int.Conference on Ultrawide-band Communications2011.Digital Object Identiﬁer10.1109/TWC.2012.040412.1106971Note that AOA does not provide better estimation accuracy than the signal propagation time based methods[10].estimation of object location if combined with high-precision timing measurement techniques[11],such as ultrawideband (UWB)signaling,which allows centimeter and even sub-millimeter accuracy,see[12],[13],and Section VII.Due to such merits,the UWB range determination is an ideal candidate for short-range object location systems and also forms the basis for the localization of sensor nodes in the IEEE802.15.4a standard[14].The algorithms based on signal propagation time can be fur-ther classiﬁed into Time of Arrival(TOA)and Time Difference of Arrival(TDOA).TOA algorithms employ the information of the absolute signal travel time from the transmitter to the target and thence to the receivers.The term“TOA”can be used in two different cases:1)there is no synchronization between transmitters and receivers and then clock bias between them exist;2)there is synchronization between transmitters and receivers and then clock bias between them does not exist. In this paper,we consider the second situation with the synchronization between the transmitter and receivers.Such synchronization can be done by cable connections between the devices,or sophisticated wireless synchronization algo-rithms[15].TDOA is employed if there is no synchronization between the transmitter and the receivers.In that case,only the receivers are synchronized.Receivers do not know the signal travel time and therefore employ the difference of signal travel times between the receivers.It is intuitive that TOA has better performance than the TDOA,since the TDOA loses information about the signal departure time[7].The TDOA/TOA positioning problems can furthermore be divided into“active”and“passive”object cases.“Active”means that the object itself is the transmitter,while“passive”means that it is not the transmitter nor receiver,but a separate (reﬂecting/scattering)object that just interacts with the signal stemming from a separate transmitter2.There are numerous papers on the TOA/TDOA location estimation for“active”objects.Regarding TDOA,the two-stage method[16]and the Approximate Maximum Likelihood Estimation[17]are shown to be able to achieve the Cramer-Rao Lower Bound(CRLB)of“active”TDOA[8].As we know,the CRLB sets the lower bound of the estimation error variance of any un-biased method.Two important TOA methods of“active”object positioning are the Least-Square Method[18]and the Approximate Maximum Likelihood Es-timation Method[17],both of which achieve the CRLB of “active”TOA.“Active”object estimation methods are used, e.g,for cellular handsets,WLAN,satellite positioning,and active RFID.2The deﬁnitions of“active”and“passive”here are different from those in radar literature.In radar literature,“passive radar”does not transmit signals and only detects transmission while“active radar”transmits signals toward targets.1536-1276/12$31.00c 2012IEEE“Passive”positioning is necessary in many practical situa-tions like crime-prevention surveillance,assets tracking,and medical patient monitoring,where the target to be localized is neither transmitter nor receiver,but a separate(reﬂect-ing/scattering)object.The TDOA positioning algorithms for “passive”objects are essentially the same as for“active”objects.For TOA,however,the synchronization creates a fundamental difference between“active”and“passive”cases. Regarding the“passive”object positioning,to the best of our knowledge,no TOA algorithms have been developed.This paper aims toﬁll this gap by proposing a TOA algorithm for passive object location estimation,which furthermore achieves the CRLB of“passive”TOA.The key contributions are:•A novel,two step estimation(TSE)method for the passive TOA based location estimation.It borrows an idea from the TDOA algorithm of[16].•CRLB for passive TOA based location estimation.When the TOA measurement error is Gaussian and small,we prove that the TSE can achieve the CRLB.Besides,it is also shown that the estimated target locations by TSE are Gaussian random variables whose covariance matrix is the inverse of the Fisher Information Matrix(FIM)related to the CRLB.We also show that in typical situations the CRLB of TOA is much lower than that of TDOA.•Experimental study of the performances of TSE.With one transmitter and three receivers equipped with UWB antennas,we perform100experimental measurements with an aluminium pole as the target.After extracting the signal travel time by high-resolution algorithms,the location of the target is evaluated by TSE.We show that the variance of estimated target location by TSE is much (33dB)lower than that by the TDOA method in[16]. The remainder of this paper is organized as follows.Section II presents the architecture of positioning system.Section III derives the TSE,followed by comparison between CRLB of TOA and TDOA algorithms in Section IV.Section V analyzes the performance of TSE.Section VI presents the simulations results.Section VII evaluates the performance of TSE based on UWB measurement.Finally Section VIII draws the conclusions.Notation:Throughout this paper,a variable with“hat”ˆ•denotes the measured/estimated values,and the“bar”¯•denotes the mean value.Bold letters denote vectors/matrices. E(•)is the expectation operator.If not particularly speciﬁed,“TOA”in this paper denotes the“TOA”for a passive object.II.A RCHITECTURE OF L OCALIZATION S YSTEMIn this section,weﬁrst discuss the challenges of localization systems,and present the focus of this paper.Then,the system model of individual localization is discussed.A.Challenges for target localizationFor easy understanding,we consider an intruder localization system using UWB signals.Note that the intruder detection can also be performed using other methods such as the Device-free Passive(DfP)approach[19]and Radio Frequency Identiﬁcation(RFID)method[20].However,both the DfP and RFID methods are based on preliminary environmental measurement information like“Radio Map Construction”[19] and“ﬁngerprints”[20].On the other hand,the TOA based approach considered in our framework does not require the preliminary efforts for obtaining environmental information. With this example,we show the challenges of target po-sitioning system:Multiple Source Separation,Indirect Path Detection and Individual Target Localization.The intruder detection system localizes,and then directs a camera to capture the photo of the targets(intruders).This localization system consists of one transmitter and several receivers.The transmitter transmits signals which are reﬂected by the targets,then,the receivers localize the targets based on the received signals.Multiple Source Separation:If there are more than one intruders,the system needs to localize each of them.With multiple targets,each receiver receives impulses from several objects.Only the information(such as TOA)extracted from impulses reﬂected by the same target should be combined for localization.Thus,the Multiple Source Separation is very important for target localization and several techniques have been proposed for this purpose.In[21],a pattern recognition scheme is used to perform the Multiple Source Separation. Video imaging and blind source separation techniques are employed for target separation in[22].Indirect Path Detection:The transmitted signals are not only reﬂected by the intruders,but also by surrounding objects,such as walls and tables.To reduce the adverse impact of non-target objects in the localization of target, the localization process consists of two steps.In the initial/ﬁrst stage,the system measures and then stores the channel impulses without the intruders.These impulses are reﬂected by non-target objects,which is referred to as reﬂectors here.The radio signal paths existing without the target are called background paths.When the intruders are present,the system performs the second measurement. To obtain the impulses related to the intruders,the system subtracts the second measurement with theﬁrst one. The remaining impulses after the subtraction can be through one of the following paths:a)transmitter-intruders-receivers,b)transmitter-reﬂectors-intruders-receivers,c) transmitter-intruders-reﬂectors-receivers,d)transmitter-reﬂectors-intruders-reﬂectors-receivers3.Theﬁrst kind of paths are called direct paths and the rest are called indirect paths.In most situations,only direct paths can be used for localization.In the literature,there are several methods proposed for indirect path identiﬁcation[23],[24]. Individual Target Localization:After the Multiple Source Separation and Indirect Path Detection,the positioning system knows the signal impulses through the direct paths for each target.Then,the system extracts the characteristics of direct paths such as TOA and AOA.Based on these characteristics, the targets areﬁnally localized.Most researches on Individual Target Localization assumes that Multiple Source Separation and Indirect Path Detection are perfectly performed such as [16],[25]and[26].Note that the three challenges sometimes 3Note that here we omit the impulses having two or more interactions with the intruder because of the resulted low signal-to-noise radio(SNR)by multiple reﬂections.Cable for synchronizationFig.1.Illustration of TOA based Location Estimation System Model.are jointly addressed,so that the target locations are estimated in one step such as the method presented in [27].In this paper,we focus on the Individual Target Local-ization,under the same framework of [16],[25]and [26],assuming that Multiple Source Separation and Indirect Path Detection are perfectly performed in prior.In addition,we only use the TOA information for localization,which achieves very high accuracy with ultra-wideband signals.The method to ex-tract TOA information using background channel cancelation is described in details in [28]and also Section VII.B.System Model of Individual LocalizationFor ease of exposition,we consider the passive object (target)location estimation problem in a two-dimensional plane as shown in Fig.1.There is a target whose location [x,y ]is to be estimated by a system with one transmitter and M receivers.Without loss of generality,let the location of the transmitter be [0,0],and the location of the i th receiver be [a i ,b i ],1≤i ≤M .The transmitter transmits an impulse;the receivers subsequently receive the signal copies reﬂected from the target and other objects.We adopt the assumption also made in [16],[17]that the target reﬂects the signal into all ing (wired)backbone connections be-tween the transmitter and receivers,or high-accuracy wireless synchronization algorithms,the transmitter and receivers are synchronized.The errors of cable synchronization are negli-gible compared with the TOA measurement errors.Thus,at the estimation center,signal travel times can be obtained by comparing the departure time at the transmitter and the arrival time at the receivers.Let the TOA from the transmitter via the target to the i th receiver be t i ,and r i =c 0t i ,where c 0is the speed of light,1≤i ≤M .Then,r i = x 2+y 2+(x −a i )2+(y −b i )2i =1,...M.(1)For future use we deﬁne r =[r 1,r 2,...,r M ].Assuming each measurement involves an error,we haver i −ˆri =e i ,1≤i ≤M,where r i is the true value,ˆr i is the measured value and e i is the measurement error.In our model,the indirect paths areignored and we assume e i to be zero mean.The estimation system tries to ﬁnd the [ˆx ,ˆy ],that best ﬁts the above equations in the sense of minimizing the error varianceΔ=E [(ˆx −x )2+(ˆy −y )2].(2)Assuming the e i are Gaussian-distributed variables with zeromean and variances σ2i ,the conditional probability functionof the observations ˆr are formulated as follows:p (ˆr |z )=Ni =11√2πσi ·exp −(ˆr i −( x 2+y 2+ (x −a i )2+(y −b i )2))22σ2i,(3)where z =[x,y ].III.TSE M ETHODIn this section,we present the two steps of TSE andsummarize them in Algorithm 1.In the ﬁrst step of TSE,we assume x ,y , x 2+y 2are independent of each other,and obtain temporary results for the target location based on this assumption.In the second step,we remove the assumption and update the estimation results.A.Step 1of TSEIn the ﬁrst step of TSE,we obtain an initial estimate of[x,y, x 2+y 2],which is performed in two stages:Stage A and Stage B.The basic idea here is to utilize the linear approximation [16][29]to simplify the problem,considering that TOA measurement errors are small with UWB signals.Let v =x 2+y 2,taking the squares of both sides of (1)leads to2a i x +2b i y −2r i v =a 2i +b 2i −r 2i .Since r i −ˆr i =e i ,it follows that−a 2i +b 2i −ˆr 2i 2+a i x +b i y −ˆr i v=e i (v −ˆr i )−e 2i 2=e i (v −ˆr i )−O (e 2i ).(4)where O (•)is the Big O Notation meaning that f (α)=O (g (α))if and only if there exits a positive real number M and a real number αsuch that|f (α)|≤M |g (α)|for all α>α0.If e i is small,we can omit the second or higher order terms O (e 2i )in Eqn (4).In the following of this paper,we do this,leaving the linear (ﬁrst order)term.Since there are M such equations,we can express them in a matrix form as followsh −S θ=Be +O (e 2)≈Be ,(5)whereh=⎡⎢⎢⎢⎢⎣−a21+b21−ˆr212−a22+b22−ˆr222...−a2M+b2M−ˆr2M2⎤⎥⎥⎥⎥⎦,S=−⎡⎢⎢⎢⎣a1b1−ˆr1a2b2−ˆr2...a Mb M−ˆr M⎤⎥⎥⎥⎦,θ=[x,y,v]T,e=[e1,e2,...,e M]T,andB=v·I−diag([r1,r2,...,r M]),(6) where O(e2)=[O(e21),O(e22),...,O(e2M)]T and diag(a) denotes the diagonal matrix with elements of vector a on its diagonal.For notational convenience,we deﬁne the error vectorϕ=h−Sθ.(7) According to(5)and(7),the mean ofϕis zero,and its covariance matrix is given byΨ=E(ϕϕT)=E(Bee T B T)+E(O(e2)e T B T)+E(Be O(e2)T)+E(O(e2)O(e2)T)≈¯BQ¯B T(8)where Q=diag[σ21,σ22,...,σ2M].Because¯B depends on the true values r,which are not obtainable,we use B(derived from the measurementsˆr)in our calculations.From(5)and the deﬁnition ofϕ,it follows thatϕis a vector of Gaussian variables;thus,the probability density function (pdf)ofϕgivenθisp(ϕ|θ)≈1(2π)M2|Ψ|12exp(−12ϕTΨ−1ϕ)=1(2π)M2|Ψ|12exp(−12(h−Sθ)TΨ−1(h−Sθ)).Then,lnp(ϕ|θ)≈−12(h−Sθ)TΨ−1(h−Sθ)+ln|Ψ|−M2ln2π(9)We assume for the moment that x,y,v are independent of each other(this clearly non-fulﬁlled assumption will be relaxed in the second step of the algorithm).Then,according to(9),the optimumθthat maximizes p(ϕ|θ)is equivalent to the one minimizingΠ=(h−Sθ)TΨ−1(h−Sθ)+ln|Ψ|. IfΨis a constant,the optimumθto minimizeΠsatisﬁes dΠdθθ=0.Taking the derivative ofΠoverθ,we havedΠdθθ=−2S TΨ−1h+2S TΨ−1Sθ.Fig.2.Illustration of estimation ofθin step1of TSE.Thus,the optimumθsatisﬁesˆθ=arg minθ{Π}=(S TΨ−1S)−1S TΨ−1h,(10)which provides[ˆx,ˆy].Note that(10)also provides the leastsquares solution for non-Gaussian errors.However,for our problem,Ψis a function ofθsince Bdepends on the(unknown)values[x,y].For this reason,themaximum-likelihood(ML)estimation method in(10)can notbe directly used.Toﬁnd the optimumθ,we perform theestimation in two stages:Stage A and Stage B.In Stage A,themissing data(Ψ)is calculated given the estimate of parameters(θ).Note thatθprovides the values of[x,y]and thus thevalue of B,therefore,Ψcan be calculated usingθby(8).In the Stage B,the parameters(θ)are updated according to(10)to maximize the likelihood function(which is equivalentto minimizingΠ).These two stages are iterated until con-vergence.Simulations in Section V show that commonly oneiteration is enough for TSE to closely approach the CRLB,which indicates that the global optimum is reached.B.Step2of TSEIn the above calculations,ˆθcontains three componentsˆx,ˆy andˆv.They were previously assumed to be independent;however,ˆx andˆy are clearly not independent ofˆv.As amatter of fact,we wish to eliminateˆv;this will be achievedby treatingˆx,ˆy,andˆv as random variables,and,knowing thelinear mapping of their squared values,the problem can besolved using the LS solution.Letˆθ=⎡⎣ˆxˆyˆv⎤⎦=⎡⎣x+n1y+n2v+n3⎤⎦(11)where n i(i=1,2,3)are the estimation errors of theﬁrststep.Obviously,the estimator(10)is an unbiased one,and themean of n i is zero.Before proceeding,we need the following Lemma.Lemma 1:By omitting the second or higher order errors,the covariance of ˆθcan be approximated as cov (ˆθ)=E (nn T )≈(¯S T Ψ−1¯S )−1.(12)where n =[n 1,n 2,n 3]T ,and Ψand ¯S(the mean value of S )use the true/mean values of x ,y,and r i .Proof:Please refer to the Appendix.Note that since the true values of x ,y,and r i are not obtain-able,we use the estimated/measured values in the calculationof cov (ˆθ).Let us now construct a vector g as followsg =ˆΘ−G Υ,(13)where ˆΘ=[ˆx 2,ˆy 2,ˆv 2]T ,Υ=[x 2,y 2]T and G =⎡⎣100111⎤⎦.Note that here ˆΘis the square of estimation result ˆθfrom the ﬁrst step containing the estimated values ˆx ,ˆy and ˆv .Υis the vector to be estimated.If ˆΘis obtained without error,g =0and the location of the target is perfectly obtained.However,the error inevitably exists and we need to estimate Υ.Recalling that v =x 2+y 2,substituting (11)into (13),and omitting the second-order terms n 21,n 22,n 23,it follows that,g =⎡⎣2xn 1+O (n 21)2yn 2+O (n 22)2vn 3+O (n 23)⎤⎦≈⎡⎣2xn 12yn 22vn 3⎤⎦.Besides,following similar procedure as that in computing(8),we haveΩ=E (gg T )≈4¯D cov (ˆθ)¯D ,(14)where ¯D =diag ([¯x ,¯y ,¯v ]).Since x ,y are not known,¯Dis calculated as ˆD using the estimated values ˆx ,ˆy from the ﬁrststep.The vector g can be approximated as a vector of Gaussian variables.Thus the maximum likelihood estimation of Υis theone minimizing (ˆΘ−G Υ)T Ω−1(ˆΘ−G Υ),expressed by ˆΥ=(G T Ω−1G )−1G T Ω−1ˆΘ.(15)The value of Ωis calculated according to (14)using the valuesof ˆx and ˆy in the ﬁrst step.Finally,the estimation of target location z is obtained byˆz =[ˆx ,ˆy ]=[±ˆΥ1,± ˆΥ2],(16)where ˆΥi is the i th item of Υ,i =1,2.To choose the correct one among the four values in (16),we can test the square error as followsχ=M i =1( ˆx 2+ˆy 2+ (ˆx −a i )2+(ˆy −b i )−ˆr i )2.(17)The value of z that minimizes χis considered as the ﬁnal estimate of the target location.In summary,the procedure of TSE is listed in Algorithm 1:Note that one should avoid placing the receivers on a line,since in this case (S T Ψ−1S )−1can become nearly singular,and solving (10)is not accurate.Algorithm 1TSE Location Estimation Method1.In the ﬁrst step,use algorithm as shown in Fig.2to obtain ˆθ,2.In the second step,use the values of ˆx and ˆy from ˆθ,generate ˆΘand D ,and calculate Ω.Then,calculate the value of ˆΥby (15),3.Among the four candidate values of ˆz =[ˆx ,ˆy ]obtained by (16),choose the one minimizing (17)as the ﬁnal estimate for target location.IV.C OMPARISON OF CRLB BETWEEN TDOA AND TOA In this section,we derive the CRLB of TOA based estima-tion algorithms and show that it is much lower (can be 30dB lower)than the CRLB of TDOA algorithms.The CRLB of “active”TOA localization has been studied in [30].The “passive”localization has been studied before under the model of multistatic radar [31],[32],[33].The difference between our model and the radar model is that in our model the localization error is a function of errors of TOA measurements,while in the radar model the localization error is a function of signal SNR and waveform.The CRLB is related to the 2×2Fisher Information Matrix (FIM)[34],J ,whose components J 11,J 12,J 21,J 22are deﬁned in (18)–(20)as follows J 11=−E (∂2ln(p (ˆr |z ))∂x 2)=ΣM i =11σ2i (x −a i (x −a i )2+(y −b i )2+xx 2+y2)2,(18)J 12=J 21=−E (∂2ln(p (ˆr |z ))∂x∂y )=ΣM i =11σ2i (x −a i (x −a i )2+(y −b i )2+x x 2+y 2)×(y −b i (x −a i )2+(y −b i )2+yx 2+y 2),(19)J 22=−E (∂2ln(p (ˆr |z ))∂y 2)=ΣM i =11σ2i (y −b i (x −a i )2+(y −b i )2+yx 2+y2)2.(20)This can be expressed asJ =U T Q −1U ,(21)where Q is deﬁned after Eqn.(8),and the entries of U in the ﬁrst and second column are{U }i,1=x ¯r i −a ix 2+y 2(x −a i )2+(y −b i )2 x 2+y 2,(22)and{U }i,2=y ¯r i −b ix 2+y 2(x −a i )2+(y −b i )2 x 2+y 2,(23)with ¯r i =(x −a i )2+(y −b i )2+ x 2+y 2.The CRLB sets the lower bound for the variance of esti-mation error of TOA algorithms,which can be expressed as [34]E [(ˆx −x )2+(ˆy −y )2]≥ J −1 1,1+J −1 2,2=CRLB T OA ,(24)where ˆx and ˆy are the estimated values of x and y ,respec-tively,and J −1 i,j is the (i,j )th element of the inverse matrix of J in (21).For the TDOA estimation,its CRLB has been derived in [16].The difference of signal travel time between several receivers are considered:(x −a i )2+(y −b i )2−(x −a 1)2+(y −b 1)2=r i −r 1=l i ,2≤i ≤M.(25)Let l =[l 2,l 3,...,l M ]T ,and t be the observa-tions/measurements of l ,then,the conditional probability density function of t is p (t |z )=1(2π)(M −1)/2|Z |12×exp(−12(t −l )T Z −1(t −l )),where Z is the correlation matrix of t ,Z =E (tt T ).Then,the FIM is expressed as [16]ˇJ=ˇU T Z −1ˇU (26)where ˇUis a M −1×2matrix deﬁned as ˇU i,1=x −a i (x −a i )2+(y −b i )2−x −a 1(x −a 1)2+(y −b 1)2,ˇUi,2=y −b i (x −a i )2+(y −b i )2−y −b 1(x −a 1)2+(y −b 1)2.The CRLB sets the lower bound for the variance of esti-mation error of TDOA algorithms,which can be expressed as [34]:E [(ˆx −x )2+(ˆy −y )2]≥ ˇJ −1 1,1+ ˇJ −1 2,2=CRLB T DOA .(27)Note that the correlation matrix Q for TOA is different from the correlation matrix Z for TDOA.Assume the variance of TOA measurement at i th (1≤i ≤M )receiver is σ2i ,it follows that:Q (i,j )=σ2i i =j,0i =j.and Z (i,j )= σ21+σ2i +1i =j,σ21i =j.As an example,we consider a scenario wherethere is a transmitter at [0,0],and four receivers at [−6,2],[6.2,1.4],[1.5,4],[2,2.3].The range of the targetlocations is 1≤x ≤10,1≤y ≤10.The ratio of CRLB of TOA over that of TDOA is plotted in Fig.3.Fig.3(a)shows the contour plot while Fig.3(b)shows the color-coded plot.It can be observed that the CRLB of TOA is always —in most cases signiﬁcantly —lower than that of TDOA.xy(a )xy0.10.20.30.40.50.60.70.80.9Fig.3.CRLB ratio of passive TOA over passive TDOA estimation:(a)contour plot;(b)pcolor plot.V.P ERFORMANCE OF TSEIn this section,we ﬁrst prove that the TSE can achieve the CRLB of TOA algorithms by showing that the estimation error variance of TSE is the same as the CRLB of TOA algorithms.In addition,we show that,for small TOA error regions,the estimated target location is approximately a Gaussian random variable whose covariance matrix is the inverse of the Fisher Information Matrix (FIM),which in turn is related to the CRLB.Similar to the reasoning in Lemma 1,we can obtain the variance of error in the estimation of Υas follows:cov (ˆΥ)≈(G T Ω−1G )−1.(28)Let ˆx =x +e x ,ˆy=y +e y ,and insert them into Υ,omitting the second order errors,we obtainˆΥ1−x 2=2xe x +O (e 2x )≈2xe x ˆΥ2−y 2=2ye y +O (e 2y)≈2ye y (29)Then,the variance of the ﬁnal estimate of target location ˆzis cov (ˆz )=E (e x e ye x e y )≈14C −1E ( Υ1−x 2Υ2−y 2Υ1−x 2Υ2−y 2 )C −1=14C −1cov (ˆΥ)C −1,(30)where C = x 00y.Substituting (14),(28),(12)and (8)into (30),we can rewrite cov (ˆz )as cov (ˆz )≈(W T Q −1W )−1(31)where W =B −1¯SD−1GC .Since we are computing an error variance,B (19),¯S(5)and D (14)are calculated using the true (mean)value of x ,y and r i .Using (19)and (1),we can rewrite B =−diag ([d 1,d 2,...,d M ]),whered i=(x−a i)2+(y−b i)2.Then B−1¯SD−1is given by B−1¯SD−1=⎡⎢⎢⎢⎢⎢⎣a1xd1b1yd1−¯r1√x2+y2d1a2xd2b2yd2−¯r2√x2+y2d2.........a Mxd Mb Myd M−¯r M√x2+y2d M⎤⎥⎥⎥⎥⎥⎦.(32)Consequently,we obtain the entries of W as{W}i,1=x¯r i−a ix2+y2(x−a i)2+(y−b i)2x2+y2,(33){W}i,2=y¯r i −b ix2+y2(x−a i)2+(y−b i)2x2+y2.(34)where{W}i,j denotes the entry at the i th row and j th column.From this we can see that W=paring(21)and (31),it followscov(ˆz)≈J−1.(35) Then,E[(ˆx−x)2+(ˆy−y)2]≈J−11,1+J−12,2.Therefore,the variance of the estimation error is the same as the CRLB.In the following,weﬁrst employ an example to show that[ˆx,ˆy]obtained by TSE are Gaussian distributed with covariance matrix J−1,and then give the explanation for this phenomenon.Let the transmitter be at[0,0],target at[0.699, 4.874]and four receivers at[-1,1],[2,1],[-31.1]and[4 0].The signal travel distance variance at four receivers are [0.1000,0.1300,0.1200,0.0950]×10−4.The two dimensional probability density function(PDF)of[ˆx,ˆy]is shown in Fig.4 (a).To verify the Gaussianity of[ˆx,ˆy],the difference between the PDF of[ˆx,ˆy]and the PDF of Gaussian distribution with mean[¯x,¯y]and covariance J−1is plotted in Fig.4(b).The Gaussianity of[ˆx,ˆy]can be explained as follows.Eqn.(35)means that the covariance of theﬁnal estimation of target location is the FIM related to CRLB.We could further study the distribution of[e x,e y].The basic idea is that by omitting the second or high order and nonlinear errors,[e x,e y]can be written as linear function of e:1)According to(29),[e x,e y]are approximately lineartransformations ofˆΥ.2)(15)means thatˆΥis approximately a linear transfor-mation ofˆΘ.Here we could omit the nonlinear errors occurred in the estimate/calculation ofΩ.3)According to(11),ˆΘ≈¯θ2+2¯θn+n2,thus,omittingthe second order error,thus,ˆΘis approximately a linear transformation of n.4)(10)and(39)mean that n is approximately a lineartransformation of e.Here we could omit the nonlinear errors accrued in the estimate of S andΨ.Thus,we could approximately write[e x,e y]as a linear trans-formation of e,thus,[e x,e y]can be approximated as Gaussian variables.Fig.4.(a):PDF of[ˆx,ˆy]by TSE(b):difference between the PDF of[ˆx,ˆy] by TSE and PDF of Gaussian distribution with mean[¯x,¯y]and covariance J−1.Fig.5.Simulation results of TSE for theﬁrst conﬁguration.VI.S IMULATION R ESULTSIn this section,weﬁrst compare the performance of TSE with that TDOA algorithm proposed in[16]and CRLBs.Then, we show the performance of TSE at high TOA measurement error scenario.For comparison,the performance of a Quasi-Newton iterative method[35]is shown.To verify our theoretical analysis,six different system con-ﬁgurations are simulated.The transmitter is at[0,0]for all six conﬁgurations,and the receiver locations and error variances are listed in Table I.Figures5,6and7show simulation results comparing the distance to the target(Conﬁguration1vs. Conﬁguration2),the receiver separation(Conﬁguration3vs. Conﬁguration4)and the number of receivers(Conﬁguration5 vs.Conﬁguration6),respectively4.In eachﬁgure,10000trails are simulated and the estimation variance of TSE estimate is compared with the CRLB of TDOA and TOA based localization schemes.For comparison,the simulation results of error variance of the TDOA method proposed in[16]are also drawn in eachﬁgure.It can be observed that1)The localization error of TSE can closely approach theCRLB of TOA based positioning algorithms.4During the simulations,only one iteration is used for the calculation of B(19).。

An Overview of Recent Progress in the Study of Distributed Multi-agent Coordination

An Overview of Recent Progress in the Study of Distributed Multi-agent CoordinationYongcan Cao,Member,IEEE,Wenwu Yu,Member,IEEE,Wei Ren,Member,IEEE,and Guanrong Chen,Fellow,IEEEAbstract—This article reviews some main results and progress in distributed multi-agent coordination,focusing on papers pub-lished in major control systems and robotics journals since 2006.Distributed coordination of multiple vehicles,including unmanned aerial vehicles,unmanned ground vehicles and un-manned underwater vehicles,has been a very active research subject studied extensively by the systems and control community. The recent results in this area are categorized into several directions,such as consensus,formation control,optimization, and estimation.After the review,a short discussion section is included to summarize the existing research and to propose several promising research directions along with some open problems that are deemed important for further investigations.Index Terms—Distributed coordination,formation control,sen-sor networks,multi-agent systemI.I NTRODUCTIONC ONTROL theory and practice may date back to thebeginning of the last century when Wright Brothers attempted theirﬁrst testﬂight in1903.Since then,control theory has gradually gained popularity,receiving more and wider attention especially during the World War II when it was developed and applied toﬁre-control systems,missile nav-igation and guidance,as well as various electronic automation devices.In the past several decades,modern control theory was further advanced due to the booming of aerospace technology based on large-scale engineering systems.During the rapid and sustained development of the modern control theory,technology for controlling a single vehicle, albeit higher-dimensional and complex,has become relatively mature and has produced many effective tools such as PID control,adaptive control,nonlinear control,intelligent control, This work was supported by the National Science Foundation under CAREER Award ECCS-1213291,the National Natural Science Foundation of China under Grant No.61104145and61120106010,the Natural Science Foundation of Jiangsu Province of China under Grant No.BK2011581,the Research Fund for the Doctoral Program of Higher Education of China under Grant No.20110092120024,the Fundamental Research Funds for the Central Universities of China,and the Hong Kong RGC under GRF Grant CityU1114/11E.The work of Yongcan Cao was supported by a National Research Council Research Associateship Award at AFRL.Y.Cao is with the Control Science Center of Excellence,Air Force Research Laboratory,Wright-Patterson AFB,OH45433,USA.W.Yu is with the Department of Mathematics,Southeast University,Nanjing210096,China and also with the School of Electrical and Computer Engineering,RMIT University,Melbourne VIC3001,Australia.W.Ren is with the Department of Electrical Engineering,University of California,Riverside,CA92521,USA.G.Chen is with the Department of Electronic Engineering,City University of Hong Kong,Hong Kong SAR,China.Copyright(c)2009IEEE.Personal use of this material is permitted. However,permission to use this material for any other purposes must be obtained from the IEEE by sending a request to pubs-permissions@.and robust control methodologies.In the past two decades in particular,control of multiple vehicles has received increas-ing demands spurred by the fact that many beneﬁts can be obtained when a single complicated vehicle is equivalently replaced by multiple yet simpler vehicles.In this endeavor, two approaches are commonly adopted for controlling multiple vehicles:a centralized approach and a distributed approach. The centralized approach is based on the assumption that a central station is available and powerful enough to control a whole group of vehicles.Essentially,the centralized ap-proach is a direct extension of the traditional single-vehicle-based control philosophy and strategy.On the contrary,the distributed approach does not require a central station for control,at the cost of becoming far more complex in structure and organization.Although both approaches are considered practical depending on the situations and conditions of the real applications,the distributed approach is believed more promising due to many inevitable physical constraints such as limited resources and energy,short wireless communication ranges,narrow bandwidths,and large sizes of vehicles to manage and control.Therefore,the focus of this overview is placed on the distributed approach.In distributed control of a group of autonomous vehicles,the main objective typically is to have the whole group of vehicles working in a cooperative fashion throughout a distributed pro-tocol.Here,cooperative refers to a close relationship among all vehicles in the group where information sharing plays a central role.The distributed approach has many advantages in achieving cooperative group performances,especially with low operational costs,less system requirements,high robustness, strong adaptivity,andﬂexible scalability,therefore has been widely recognized and appreciated.The study of distributed control of multiple vehicles was perhapsﬁrst motivated by the work in distributed comput-ing[1],management science[2],and statistical physics[3]. In the control systems society,some pioneering works are generally referred to[4],[5],where an asynchronous agree-ment problem was studied for distributed decision-making problems.Thereafter,some consensus algorithms were studied under various information-ﬂow constraints[6]–[10].There are several journal special issues on the related topics published af-ter2006,including the IEEE Transactions on Control Systems Technology(vol.15,no.4,2007),Proceedings of the IEEE (vol.94,no.4,2007),ASME Journal of Dynamic Systems, Measurement,and Control(vol.129,no.5,2007),SIAM Journal of Control and Optimization(vol.48,no.1,2009),and International Journal of Robust and Nonlinear Control(vol.21,no.12,2011).In addition,there are some recent reviewsand progress reports given in the surveys[11]–[15]and thebooks[16]–[23],among others.This article reviews some main results and recent progressin distributed multi-agent coordination,published in majorcontrol systems and robotics journals since2006.Due to space limitations,we refer the readers to[24]for a more completeversion of the same overview.For results before2006,thereaders are referred to[11]–[14].Speciﬁcally,this article reviews the recent research resultsin the following directions,which are not independent but actually may have overlapping to some extent:1.Consensus and the like(synchronization,rendezvous).Consensus refers to the group behavior that all theagents asymptotically reach a certain common agreementthrough a local distributed protocol,with or without predeﬁned common speed and orientation.2.Distributed formation and the like(ﬂocking).Distributedformation refers to the group behavior that all the agents form a pre-designed geometrical conﬁguration throughlocal interactions with or without a common reference.3.Distributed optimization.This refers to algorithmic devel-opments for the analysis and optimization of large-scaledistributed systems.4.Distributed estimation and control.This refers to dis-tributed control design based on local estimation aboutthe needed global information.The rest of this article is organized as follows.In Section II,basic notations of graph theory and stochastic matrices are introduced.Sections III,IV,V,and VI describe the recentresearch results and progress in consensus,formation control, optimization,and estimation.Finally,the article is concludedby a short section of discussions with future perspectives.II.P RELIMINARIESA.Graph TheoryFor a system of n connected agents,its network topology can be modeled as a directed graph denoted by G=(V,W),where V={v1,v2,···,v n}and W⊆V×V are,respectively, the set of agents and the set of edges which directionallyconnect the agents together.Speciﬁcally,the directed edgedenoted by an ordered pair(v i,v j)means that agent j can access the state information of agent i.Accordingly,agent i is a neighbor of agent j.A directed path is a sequence of directed edges in the form of(v1,v2),(v2,v3),···,with all v i∈V.A directed graph has a directed spanning tree if there exists at least one agent that has a directed path to every other agent.The union of a set of directed graphs with the same setof agents,{G i1,···,G im},is a directed graph with the sameset of agents and its set of edges is given by the union of the edge sets of all the directed graphs G ij,j=1,···,m.A complete directed graph is a directed graph in which each pair of distinct agents is bidirectionally connected by an edge,thus there is a directed path from any agent to any other agent in the network.Two matrices are used to represent the network topology: the adjacency matrix A=[a ij]∈R n×n with a ij>0if (v j,v i)∈W and a ij=0otherwise,and the Laplacian matrix L=[ℓij]∈R n×n withℓii= n j=1a ij andℓij=−a ij,i=j, which is generally asymmetric for directed graphs.B.Stochastic MatricesA nonnegative square matrix is called(row)stochastic matrix if its every row is summed up to one.The product of two stochastic matrices is still a stochastic matrix.A row stochastic matrix P∈R n×n is called indecomposable and aperiodic if lim k→∞P k=1y T for some y∈R n[25],where 1is a vector with all elements being1.III.C ONSENSUSConsider a group of n agents,each with single-integrator kinematics described by˙x i(t)=u i(t),i=1,···,n,(1) where x i(t)and u i(t)are,respectively,the state and the control input of the i th agent.A typical consensus control algorithm is designed asu i(t)=nj=1a ij(t)[x j(t)−x i(t)],(2)where a ij(t)is the(i,j)th entry of the corresponding ad-jacency matrix at time t.The main idea behind(2)is that each agent moves towards the weighted average of the states of its neighbors.Given the switching network pattern due to the continuous motions of the dynamic agents,coupling coefﬁcients a ij(t)in(2),hence the graph topologies,are generally time-varying.It is shown in[9],[10]that consensus is achieved if the underlying directed graph has a directed spanning tree in some jointly fashion in terms of a union of its time-varying graph topologies.The idea behind consensus serves as a fundamental principle for the design of distributed multi-agent coordination algo-rithms.Therefore,investigating consensus has been a main research direction in the study of distributed multi-agent co-ordination.To bridge the gap between the study of consensus algorithms and many physical properties inherited in practical systems,it is necessary and meaningful to study consensus by considering many practical factors,such as actuation,control, communication,computation,and vehicle dynamics,which characterize some important features of practical systems.This is the main motivation to study consensus.In the following part of the section,an overview of the research progress in the study of consensus is given,regarding stochastic network topologies and dynamics,complex dynamical systems,delay effects,and quantization,mainly after2006.Several milestone results prior to2006can be found in[2],[4]–[6],[8]–[10], [26].A.Stochastic Network Topologies and DynamicsIn multi-agent systems,the network topology among all vehicles plays a crucial role in determining consensus.The objective here is to explicitly identify necessary and/or sufﬁ-cient conditions on the network topology such that consensus can be achieved under properly designed algorithms.It is often reasonable to consider the case when the network topology is deterministic under ideal communication chan-nels.Accordingly,main research on the consensus problem was conducted under a deterministicﬁxed/switching network topology.That is,the adjacency matrix A(t)is deterministic. Some other times,when considering random communication failures,random packet drops,and communication channel instabilities inherited in physical communication channels,it is necessary and important to study consensus problem in the stochastic setting where a network topology evolves according to some random distributions.That is,the adjacency matrix A(t)is stochastically evolving.In the deterministic setting,consensus is said to be achieved if all agents eventually reach agreement on a common state. In the stochastic setting,consensus is said to be achieved almost surely(respectively,in mean-square or in probability)if all agents reach agreement on a common state almost surely (respectively,in mean-square or with probability one).Note that the problem studied in the stochastic setting is slightly different from that studied in the deterministic setting due to the different assumptions in terms of the network topology. Consensus over a stochastic network topology was perhaps ﬁrst studied in[27],where some sufﬁcient conditions on the network topology were given to guarantee consensus with probability one for systems with single-integrator kinemat-ics(1),where the rate of convergence was also studied.Further results for consensus under a stochastic network topology were reported in[28]–[30],where research effort was conducted for systems with single-integrator kinematics[28],[29]or double-integrator dynamics[30].Consensus for single-integrator kine-matics under stochastic network topology has been exten-sively studied in particular,where some general conditions for almost-surely consensus was derived[29].Loosely speaking, almost-surely consensus for single-integrator kinematics can be achieved,i.e.,x i(t)−x j(t)→0almost surely,if and only if the expectation of the network topology,namely,the network topology associated with expectation E[A(t)],has a directed spanning tree.It is worth noting that the conditions are analogous to that in[9],[10],but in the stochastic setting. In view of the special structure of the closed-loop systems concerning consensus for single-integrator kinematics,basic properties of the stochastic matrices play a crucial role in the convergence analysis of the associated control algorithms. Consensus for double-integrator dynamics was studied in[30], where the switching network topology is assumed to be driven by a Bernoulli process,and it was shown that consensus can be achieved if the union of all the graphs has a directed spanning tree.Apparently,the requirement on the network topology for double-integrator dynamics is a special case of that for single-integrator kinematics due to the difference nature of theﬁnal states(constantﬁnal states for single-integrator kinematics and possible dynamicﬁnal states for double-integrator dynamics) caused by the substantial dynamical difference.It is still an open question as if some general conditions(corresponding to some speciﬁc algorithms)can be found for consensus with double-integrator dynamics.In addition to analyzing the conditions on the network topology such that consensus can be achieved,a special type of consensus algorithm,the so-called gossip algorithm[31],[32], has been used to achieve consensus in the stochastic setting. The gossip algorithm can always guarantee consensus almost surely if the available pairwise communication channels satisfy certain conditions(such as a connected graph).The way of network topology switching does not play any role in the consideration of consensus.The current study on consensus over stochastic network topologies has shown some interesting results regarding:(1) consensus algorithm design for various multi-agent systems,(2)conditions of the network topologies on consensus,and(3)effects of the stochastic network topologies on the con-vergence rate.Future research on this topic includes,but not limited to,the following two directions:(1)when the network topology itself is stochastic,how to determine the probability of reaching consensus almost surely?(2)compared with the deterministic network topology,what are the advantages and disadvantages of the stochastic network topology,regarding such as robustness and convergence rate?As is well known,disturbances and uncertainties often exist in networked systems,for example,channel noise,commu-nication noise,uncertainties in network parameters,etc.In addition to the stochastic network topologies discussed above, the effect of stochastic disturbances[33],[34]and uncertain-ties[35]on the consensus problem also needs investigation. Study has been mainly devoted to analyzing the performance of consensus algorithms subject to disturbances and to present-ing conditions on the uncertainties such that consensus can be achieved.In addition,another interesting direction in dealing with disturbances and uncertainties is to design distributed localﬁltering algorithms so as to save energy and improve computational efﬁciency.Distributed localﬁltering algorithms play an important role and are more effective than traditional centralizedﬁltering algorithms for multi-agent systems.For example,in[36]–[38]some distributed Kalmanﬁlters are designed to implement data fusion.In[39],by analyzing consensus and pinning control in synchronization of complex networks,distributed consensusﬁltering in sensor networks is addressed.Recently,Kalmanﬁltering over a packet-dropping network is designed through a probabilistic approach[40]. Today,it remains a challenging problem to incorporate both dynamics of consensus and probabilistic(Kalman)ﬁltering into a uniﬁed framework.plex Dynamical SystemsSince consensus is concerned with the behavior of a group of vehicles,it is natural to consider the system dynamics for practical vehicles in the study of the consensus problem. Although the study of consensus under various system dynam-ics is due to the existence of complex dynamics in practical systems,it is also interesting to observe that system dynamics play an important role in determining theﬁnal consensus state.For instance,the well-studied consensus of multi-agent systems with single-integrator kinematics often converges to a constantﬁnal value instead.However,consensus for double-integrator dynamics might admit a dynamicﬁnal value(i.e.,a time function).These important issues motivate the study of consensus under various system dynamics.As a direct extension of the study of the consensus prob-lem for systems with simple dynamics,for example,with single-integrator kinematics or double-integrator dynamics, consensus with general linear dynamics was also studied recently[41]–[43],where research is mainly devoted toﬁnding feedback control laws such that consensus(in terms of the output states)can be achieved for general linear systems˙x i=Ax i+Bu i,y i=Cx i,(3) where A,B,and C are constant matrices with compatible sizes.Apparently,the well-studied single-integrator kinematics and double-integrator dynamics are special cases of(3)for properly choosing A,B,and C.As a further extension,consensus for complex systems has also been extensively studied.Here,the term consensus for complex systems is used for the study of consensus problem when the system dynamics are nonlinear[44]–[48]or with nonlinear consensus algorithms[49],[50].Examples of the nonlinear system dynamics include:•Nonlinear oscillators[45].The dynamics are often as-sumed to be governed by the Kuramoto equation˙θi=ωi+Kstability.A well-studied consensus algorithm for(1)is given in(2),where it is now assumed that time delay exists.Two types of time delays,communication delay and input delay, have been considered in the munication delay accounts for the time for transmitting information from origin to destination.More precisely,if it takes time T ij for agent i to receive information from agent j,the closed-loop system of(1)using(2)under aﬁxed network topology becomes˙x i(t)=nj=1a ij(t)[x j(t−T ij)−x i(t)].(7)An interpretation of(7)is that at time t,agent i receives information from agent j and uses data x j(t−T ij)instead of x j(t)due to the time delay.Note that agent i can get its own information instantly,therefore,input delay can be considered as the summation of computation time and execution time. More precisely,if the input delay for agent i is given by T p i, then the closed-loop system of(1)using(2)becomes˙x i(t)=nj=1a ij(t)[x j(t−T p i)−x i(t−T p i)].(8)Clearly,(7)refers to the case when only communication delay is considered while(8)refers to the case when only input delay is considered.It should be emphasized that both communication delay and input delay might be time-varying and they might co-exist at the same time.In addition to time delay,it is also important to consider packet drops in exchanging state information.Fortunately, consensus with packet drops can be considered as a special case of consensus with time delay,because re-sending packets after they were dropped can be easily done but just having time delay in the data transmission channels.Thus,the main problem involved in consensus with time delay is to study the effects of time delay on the convergence and performance of consensus,referred to as consensusabil-ity[52].Because time delay might affect the system stability,it is important to study under what conditions consensus can still be guaranteed even if time delay exists.In other words,can oneﬁnd conditions on the time delay such that consensus can be achieved?For this purpose,the effect of time delay on the consensusability of(1)using(2)was investigated.When there exists only(constant)input delay,a sufﬁcient condition on the time delay to guarantee consensus under aﬁxed undirected interaction graph is presented in[8].Speciﬁcally,an upper bound for the time delay is derived under which consensus can be achieved.This is a well-expected result because time delay normally degrades the system performance gradually but will not destroy the system stability unless the time delay is above a certain threshold.Further studies can be found in, e.g.,[53],[54],which demonstrate that for(1)using(2),the communication delay does not affect the consensusability but the input delay does.In a similar manner,consensus with time delay was studied for systems with different dynamics, where the dynamics(1)are replaced by other more complex ones,such as double-integrator dynamics[55],[56],complex networks[57],[58],rigid bodies[59],[60],and general nonlinear dynamics[61].In summary,the existing study of consensus with time delay mainly focuses on analyzing the stability of consensus algo-rithms with time delay for various types of system dynamics, including linear and nonlinear dynamics.Generally speaking, consensus with time delay for systems with nonlinear dynam-ics is more challenging.For most consensus algorithms with time delays,the main research question is to determine an upper bound of the time delay under which time delay does not affect the consensusability.For communication delay,it is possible to achieve consensus under a relatively large time delay threshold.A notable phenomenon in this case is that theﬁnal consensus state is constant.Considering both linear and nonlinear system dynamics in consensus,the main tools for stability analysis of the closed-loop systems include matrix theory[53],Lyapunov functions[57],frequency-domain ap-proach[54],passivity[58],and the contraction principle[62]. Although consensus with time delay has been studied extensively,it is often assumed that time delay is either constant or random.However,time delay itself might obey its own dynamics,which possibly depend on the communication distance,total computation load and computation capability, etc.Therefore,it is more suitable to represent the time delay as another system variable to be considered in the study of the consensus problem.In addition,it is also important to consider time delay and other physical constraints simultaneously in the study of the consensus problem.D.QuantizationQuantized consensus has been studied recently with motiva-tion from digital signal processing.Here,quantized consensus refers to consensus when the measurements are digital rather than analog therefore the information received by each agent is not continuous and might have been truncated due to digital ﬁnite precision constraints.Roughly speaking,for an analog signal s,a typical quantizer with an accuracy parameterδ, also referred to as quantization step size,is described by Q(s)=q(s,δ),where Q(s)is the quantized signal and q(·,·) is the associated quantization function.For instance[63],a quantizer rounding a signal s to its nearest integer can be expressed as Q(s)=n,if s∈[(n−1/2)δ,(n+1/2)δ],n∈Z, where Z denotes the integer set.Note that the types of quantizers might be different for different systems,hence Q(s) may differ for different systems.Due to the truncation of the signals received,consensus is now considered achieved if the maximal state difference is not larger than the accuracy level associated with the whole system.A notable feature for consensus with quantization is that the time to reach consensus is usuallyﬁnite.That is,it often takes aﬁnite period of time for all agents’states to converge to an accuracy interval.Accordingly,the main research is to investigate the convergence time associated with the proposed consensus algorithm.Quantized consensus was probablyﬁrst studied in[63], where a quantized gossip algorithm was proposed and its convergence was analyzed.In particular,the bound of theconvergence time for a complete graph was shown to be poly-nomial in the network size.In[64],coding/decoding strate-gies were introduced to the quantized consensus algorithms, where it was shown that the convergence rate depends on the accuracy of the quantization but not the coding/decoding schemes.In[65],quantized consensus was studied via the gossip algorithm,with both lower and upper bounds of the expected convergence time in the worst case derived in terms of the principle submatrices of the Laplacian matrix.Further results regarding quantized consensus were reported in[66]–[68],where the main research was also on the convergence time for various proposed quantized consensus algorithms as well as the quantization effects on the convergence time.It is intuitively reasonable that the convergence time depends on both the quantization level and the network topology.It is then natural to ask if and how the quantization methods affect the convergence time.This is an important measure of the robustness of a quantized consensus algorithm(with respect to the quantization method).Note that it is interesting but also more challenging to study consensus for general linear/nonlinear systems with quantiza-tion.Because the difference between the truncated signal and the original signal is bounded,consensus with quantization can be considered as a special case of one without quantization when there exist bounded disturbances.Therefore,if consensus can be achieved for a group of vehicles in the absence of quantization,it might be intuitively correct to say that the differences among the states of all vehicles will be bounded if the quantization precision is small enough.However,it is still an open question to rigorously describe the quantization effects on consensus with general linear/nonlinear systems.E.RemarksIn summary,the existing research on the consensus problem has covered a number of physical properties for practical systems and control performance analysis.However,the study of the consensus problem covering multiple physical properties and/or control performance analysis has been largely ignored. In other words,two or more problems discussed in the above subsections might need to be taken into consideration simul-taneously when studying the consensus problem.In addition, consensus algorithms normally guarantee the agreement of a team of agents on some common states without taking group formation into consideration.To reﬂect many practical applications where a group of agents are normally required to form some preferred geometric structure,it is desirable to consider a task-oriented formation control problem for a group of mobile agents,which motivates the study of formation control presented in the next section.IV.F ORMATION C ONTROLCompared with the consensus problem where theﬁnal states of all agents typically reach a singleton,theﬁnal states of all agents can be more diversiﬁed under the formation control scenario.Indeed,formation control is more desirable in many practical applications such as formationﬂying,co-operative transportation,sensor networks,as well as combat intelligence,surveillance,and reconnaissance.In addition,theperformance of a team of agents working cooperatively oftenexceeds the simple integration of the performances of all individual agents.For its broad applications and advantages,formation control has been a very active research subject inthe control systems community,where a certain geometric pattern is aimed to form with or without a group reference.More precisely,the main objective of formation control is to coordinate a group of agents such that they can achievesome desired formation so that some tasks can beﬁnished bythe collaboration of the agents.Generally speaking,formation control can be categorized according to the group reference.Formation control without a group reference,called formationproducing,refers to the algorithm design for a group of agents to reach some pre-desired geometric pattern in the absenceof a group reference,which can also be considered as the control objective.Formation control with a group reference,called formation tracking,refers to the same task but followingthe predesignated group reference.Due to the existence of the group reference,formation tracking is usually much morechallenging than formation producing and control algorithmsfor the latter might not be useful for the former.As of today, there are still many open questions in solving the formationtracking problem.The following part of the section reviews and discussesrecent research results and progress in formation control, including formation producing and formation tracking,mainlyaccomplished after2006.Several milestone results prior to 2006can be found in[69]–[71].A.Formation ProducingThe existing work in formation control aims at analyzingthe formation behavior under certain control laws,along with stability analysis.1)Matrix Theory Approach:Due to the nature of multi-agent systems,matrix theory has been frequently used in thestability analysis of their distributed coordination.Note that consensus input to each agent(see e.g.,(2))isessentially a weighted average of the differences between the states of the agent’s neighbors and its own.As an extensionof the consensus algorithms,some coupling matrices wereintroduced here to offset the corresponding control inputs by some angles[72],[73].For example,given(1),the controlinput(2)is revised as u i(t)= n j=1a ij(t)C[x j(t)−x i(t)], where C is a coupling matrix with compatible size.If x i∈R3, then C can be viewed as the3-D rotational matrix.The mainidea behind the revised algorithm is that the original controlinput for reaching consensus is now rotated by some angles. The closed-loop system can be expressed in a vector form, whose stability can be determined by studying the distribution of the eigenvalues of a certain transfer matrix.Main research work was conducted in[72],[73]to analyze the collective motions for systems with single-integrator kinematics and double-integrator dynamics,where the network topology,the damping gain,and C were shown to affect the collective motions.Analogously,the collective motions for a team of nonlinear self-propelling agents were shown to be affected by。

遥感与地理信息系统方面的好的期刊

遥感与地理信息系统方面的好杂志国内的期刊：1）遥感学报（98年前《环境遥感》杂志，国内比较好的遥感专业杂志，主编是原遥感所所长、现国家科技部部长徐冠华院士，遥感文章比较多，象国内比较牛的遥感理论研究的大牛复旦大学的金亚秋教授和北京师范大学的新当选的院士李小文教授经常有文章发表；基于遥感和GIS资源环境等应用的文章也比较好，主要是中科院地理所和遥感所的；还有就是图像处理的算法研究或新型的遥感方法如雷达干涉测量、高光谱方面的研究，主要由武汉大学测绘遥感信息工程国家重点实验室(L)和中科院遥感所的文章。

（2）测绘学报（侧重测量基础理论的研究，但经常有非常好的综述型的文章，上面的测绘学博士论文摘要是非常好，还有主编陈俊勇院士治学非常严谨，一般的假冒伪劣文章很难找到市场，该刊宁缺勿滥,2001年仍然是季刊,文章少,但很精。

不过该刊刊登的文章比较偏重大地测量（GPS），GIS的文章相比比较少）。

（3）武测学报（2001年改名《武汉大学学报》信息科学版）本杂志是原武汉测绘科技大学学报，主编是中国科学院和中国工程院双院士李德仁教授,很多具有创新性和理论性的测绘研究成果都在该刊发表，展示中国测绘学术研究的最高水平,引导测绘理论研究的方向。

我认为上面的博士论文摘要比较好，真正体现了我国3S技术的研究动向和学术水准。

本刊出版内容包括综述与瞻望、学术论文和研究报告、本领域重大科技新闻等,涉及测绘学研究的主要方面,尤其是摄影测量与遥感、大地测量与地球重力场、工程测量、地图制图、地球动力学、全球定位系统(GPS)、地理信息系统(GIS)、图形图像处理等。

该刊现同时有英文版，名为GEO-SPATIAL INFORMATION SCIENCE，是中文版的精华版，万方科技期刊网上可以下载全文。

（4）中国图象图形学报1996年创刊，由中国图象图形学会、中科院遥感所、中科院计算所共同主办，主编是科技部部长徐冠华院士。

2001年起《中国图象图形学报》分A、B版。

Target Detection and Localization Using MIMO

IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 54, NO. 10, OCTOBER 20063873Target Detection and Localization Using MIMO Radars and SonarsIlya Bekkerman and Joseph Tabrikian, Senior Member, IEEEAbstract—In this paper, we propose a new space-time coding conﬁguration for target detection and localization by radar or sonar systems. In common active array systems, the transmitted signal is usually coherent between the different elements of the array. This conﬁguration does not allow array processing in the transmit mode. However, space-time coding of the transmitted signals allows to digitally steer the beam pattern in the transmit in addition to the received signal. The ability to steer the transmitted beam pattern, helps to avoid beam shape loss. We show that the conﬁguration with spatially orthogonal signal transmission is equivalent to additional virtual sensors which extend the array aperture with virtual spatial tapering. These virtual sensors can be used to form narrower beams with lower sidelobes and, therefore, provide higher performance in target detection, angular estimation accuracy, and angular resolution. The generalized likelihood ratio test for target detection and the maximum likelihood and Cramér–Rao bound for target direction estimation are derived for an arbitrary signal coherence matrix. It is shown that the optimal performance is achieved for orthogonal transmitted signals. Target detection and localization performances are evaluated and studied theoretically and via simulations. Index Terms—Cramér–Rao bound (CRB), generalized likelihood ratio test (GLRT), maximum likelihood, MIMO radars, MIMO sonars, orthogonal signal transmission, space–time coding, transmit beamforming, virtual sensors.I. INTRODUCTIONACTIVE target detection and localization systems, such as radars or active sonars, usually transmit a directional beam, and the target echo signal is processed in the receive mode. In the last two decades, array processing of the received signal has been intensively investigated (see, for example, [1]). However, this conﬁguration does not allow array processing in the transmit mode. In fact, the transmission is usually performed using the phased array technique or other beam steering methods. Array processing in both transmit and receive modes is possible when the transmitted signals are spatially coded, i.e., spatially orthogonal signals. This paper addresses the problem of target detection and localization by active array using spatially coded signals. Transmission of orthogonal signals from an array is commonly used in communication systems [2]. Passive localization of orthogonal signals with known waveforms was investigatedManuscript received June 19, 2004; revised October 16, 2005. The associate editor coordinating the review of this manuscript and approving it for publication was Prof. Mats Viberg. The authors are with the Department of Electrical and Computer Engineering, Ben-Gurion University of the Negev, Beer-Sheva 84105, Israel (e-mail: ilyabek@ee.bgu.ac.il; joseph@ee.bgu.ac.il). Digital Object Identiﬁer 10.1109/TSP.2006.879267in [3]. In [4], it is shown that the conventional conﬁguration of one transmitter and two receivers and an alternative conﬁguration of two transmitters and one receiver are equivalent in terms of Cramér–Rao bound (CRB) on bearing estimation. This conﬁguration requires radiating two orthogonal signals from two transmitters. The potential advantage of this conﬁguration over the conventional one is in applications where the receiving elements are to be placed on a platform of limited size. The results in [4] were extended in [5], in which three possible combinations of four transmitters/receivers were investigated: 1) one transmitter and three receivers, 2) two transmitters and two receivers, and 3) three transmitters and one receiver. It was found that these conﬁgurations have identical performance in terms of angle estimation accuracy, where the transmitting signals are orthogonal. In [6], spatio–temporal coding for an antenna array was introduced, and it was shown that a single receiver is sufﬁcient for digital beamforming. Fishler et al. [7] also investigated the problem of orthogonal signal transmission for multiple-input multiple-output (MIMO) radar. They assumed a multistatic radar in which the spacing between the elements of the array is very large, and that the transmitter and the receiver of the radar are separated such that they experience an angular spread. In [8], a novel conﬁguration for array processing using spacetime coding of the transmitted signal was presented. This conﬁguration does not assume a multistatic radar, with an arbitrary signal coherence across the radar elements. In this paper, we analyze the properties of the space-time coding conﬁguration for target detection and localization. In particular, the advantages of the proposed conﬁguration are analytically demonstrated and compared to the conventional coherent signal transmission case. The main advantages of this new conﬁguration are as follows: • digital beamforming of the transmitted beams in addition to the received beams, therefore avoiding beam shape loss in cases when the target is not in the center of the beam; • extension of the array aperture by virtual sensors, therefore obtaining narrower beams; • virtual spatial tapering of the extended array aperture, therefore obtaining lower side lobes; • improving the angular resolution by using the information in the transmit and the receive modes; • increasing the upper limit on the number of targets which can be detected and localized by the array (this is attributed to the virtual sensors); • decreasing the spatial transmitted peak power density. This paper is organized as follows. The spatially coded signal model is presented in Section II. In Section III, the sufﬁcient statistic (SS) for detection and estimation algorithms are derived. The model’s properties are analyzed in Section IV. The1053-587X/$20.00 © 2006 IEEEAuthorized licensed use limited to: NORTHWESTERN POLYTECHNIC UNIVERSITY. Downloaded on January 7, 2010 at 03:53 from IEEE Xplore. Restrictions apply.3874IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 54, NO. 10, OCTOBER 2006wherestands for the complex amplitude of the received signal, is the additive noise at the th element, and describes the total phase delay of the signal, transmitted by the th element and received by the th element, is the carrier frequency. The total delay from the th where transmitting element to the th receiving element for the farﬁeld case can be written as (3)and stands for the signal where wavelength. denote the response of the th element. Then, the Let th element of the array response can be decomposed asFig. 1. Array conﬁguration.(4) maximum likelihood (ML) estimator, the CRB for target localization, and the generalized likelihood ratio test (GLRT) for target detection are derived in Section V. The proposed concept is tested via a few examples and simulations, which appear in Section VI. The main results of this paper are discussed and concluded in Section VII. II. SPATIALLY CODED SIGNAL MODEL Consider an element antenna array transmitting narrow-band signals. The samples of baseband equivalent with coherence signals are denoted by the vectors matrix Note that the elements of depend on through all possible combinations of delays in transmit and receive modes. In fact, is the array response for transmit from the th element and receive by the th element. Hence, the array response matrix can be deﬁned as (5) in which is the transmitted or the received array response vector. In matrix notation, (2) can be written as (6) where , , and are vectors of the received signal, the transmitted signal, and the additive noise, respectively. In sensors with range or Doppler estimation capability, the model should include also the target range and Doppler. Typically, in these sensors, the ML estimator of the target direction, range, and Doppler is implemented by processing the receiving channels over time, obtaining multichannel measurements for each considered range-Doppler bin. The above model refers to a single range-Doppler bin. Speciﬁcally, the target detection and localization algorithms presented in this paper should be applied for each range-Doppler bin individually. In the above model, multiple targets can be allowed if they do not share the same range-Doppler bin. In the following, this model is extended to allow multiple targets in the considered range-Doppler bin. In the case of targets in the given range-Doppler bin, (6) is modiﬁed to (7) Let denote the vector of unknown parameters, which includes the directions of arrival (DOAs) and the complex and amplitudes of the targets, where , . The vectors and are considered as deterministic unknown.. . .. . ..... . .(1)where represents the time index, is the complex correlation coefﬁcient between the th and th signals, and denotes the Hermitian operation. The phases of control the transmitted beam direction of the coherent component. In the and the case of orthonormal transmitted signals , i.e., omnidicoherence matrix is an identity matrix: rectional transmission. In common radar systems, coherent sigis nals are transmitted by the array and therefore the rank of equal to one.1 Let denote the location of the th element of the array with denoting the transposition operation (see Fig. 1). In the presence of a single target at direction in a multipath-free environment, the received signal by the th element is given by(2)1The different elements transmit the same signal with phase shifts for beam steering.Authorized licensed use limited to: NORTHWESTERN POLYTECHNIC UNIVERSITY. Downloaded on January 7, 2010 at 03:53 from IEEE Xplore. Restrictions apply.BEKKERMAN AND TABRIKIAN: TARGET DETECTION AND LOCALIZATION USING MIMO RADARS AND SONARS3875The noise vectors are assumed to be independent, zero-mean complex Gaussian with known covariance . With no loss of generality, we can assume that matrix , where is an identity matrix of size . If this assumption is not satisﬁed, the model in (7) can be prewhitened. In many practical scenarios, in the presence of clutter or jammer, the noise covariance matrix is unknown. In these cases, the noise can be treated as a composition of additional interference sources. This problem can be solved by either multiple target localization techniques, as is modeled below, or by adaptive methods such as the Capon’s beamformer. III. SUFFICIENT STATISTIC FOR DOA ESTIMATION According to the assumptions stated in the previous section, the measurement vectors are independent complex Gaussian , where vectors with denotes the complex Gaussian distribution. The log-likelihood function for estimating from the data is derived in Appendix A and given byFig. 2. Sufﬁcient statistic extraction.The independent sufﬁcient statistic vector can be obtained as (12) The conﬁguration for obtaining the sufﬁcient statistic from the data is described in Fig. 2. Actually, the sufﬁcient statistic can be obtained by a matched ﬁlter: temporal matching the mea. surement vectors to different signal subspace components Insertion of (7) and (11) into (12) yields(8) is the where gate operation, and th column of , denotes the conjuis the th sufﬁcient statistic, deﬁned as (9) which is obtained by matching the observed data to the th . Moreover, the sufﬁcient statistic matrix can signal be deﬁned as (14) (10) Finally, (14) can be written in the form It can be shown that for nonorthogonal signals, the sufﬁcient statistic are statistically dependent. For simplicity of the algorithms, we are interested in independent sufﬁ, which can be obtained as follows. We cient statistic is nonsingular, and then will ﬁrst assume that the matrix we will refer to the general case. The matrix from (1) can be decomposed using singular value decomposition (SVD) as , where and are the matrices of eigenvec, respectively. Accordingly, is a tors and eigenvalues of linear transformation of a vector of independent signals deﬁned as (11) (15) (16) where is the equivalent array response of size at the direction and . The subscript denotes dependence of the new steering vector on the signal correlation coefﬁcients and , which is zero-mean complex Gaussian with . Actually, (15) states an equivalent model to (7).(13) By recalling the deﬁnition of from (1) and using its SVD decomposition form, (13) can be rewritten asAuthorized licensed use limited to: NORTHWESTERN POLYTECHNIC UNIVERSITY. Downloaded on January 7, 2010 at 03:53 from IEEE Xplore. Restrictions apply.3876IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 54, NO. 10, OCTOBER 2006If the matrix is singular, then part of its eigenvalues are equal to zero. In such a case, the matrix is changed to , and the above procedure can be repeated with this modiﬁed eigenvalues matrix. The ﬁnal result in this case is independent of , and therefore the “lim” can be . dropped. Accordingly, (15) and (16) hold also for singular IV. MODEL’S PROPERTIES In this section, the properties of the equivalent model of (15) are investigated and illustrated. A. Virtual Aperture Extension Recalling (16), the equivalent array response can be calculated for coherent and orthogonal signals. In the coherent has a single nonzero eigensignal case, the matrix and the equivalent value. Therefore, steering vector becomes (17) where represents a vector of size , with zero elements. By substitution of (5) into (17), one obtainsFig. 3. Array aperture for coherent signals M = 3, L = 1.(18) The equivalent array response is given by the steering vector multiplied by the gain achieved in in the receive mode the transmit mode , where is the weighting vector in the transmit mode. This gain is decreased due to the beam shape loss when the target is not located in the center of the transmit beam. In the partic, with denoting a column vector of size ular case of whose elements are equal to one, the transmit beam is directed , that is, . to the array broadside and , The equivalent array response in this case, denoted by can be written as (19) In the case of orthonormal signals, and are equal to one. Note the eigenvalues of the matrix that in this case, the matrix of eigenvectors is not unique. A simple choice for is . ThereforeFig. 4. Array aperture for orthogonal signals M = 3, L = 1; -points: actual sensors, o-points: virtual sensors. Two sensors are located at points B , C , E .2(20) The equivalent steering vector for orthonormal signals is the and product of the steering vector in the receive mode . The equivathe steering vector in the transmit mode lent steering vector for noncoherent signals includes all the elements of , which represent all the possible transmit–receive combinations. The th element of this matrix . Hence, the array is response consists of virtual sensors located at the combinations for . Consequently, the array of aperture is virtually extended. This virtual aperture extensionresults in narrower beams and therefore higher angular resolution and better detection performance. Moreover, some of the virtual sensor locations are identical, which can be interpreted as spatial tapering and results in lower sidelobes. In order to illustrate these advantages, we examine an ex, which are located ample with three array elements at vertexes of an equilateral triangle (see -points in Fig. 3) and one target . In Figs. 3 and 4, the equivalent array structure for coherent and orthogonal signals is presented, respectively. As mentioned above, the equivalent array for orthogonal signals includes all the transmit–receive combinations of the el, which are given by ements of . This is equivalent to an extended array, whose elements are located at . Therefore the equivalent array consists of virtual sensors (marked by -points in Fig. 4) in addition to the actual sensors. According to (3), the total number of delays is nine, representing nine different virtual sensors. In this conﬁguration, we obtain three sensors at the actual sensor locations (points A, B, C), three virtual sensors at new locations (points D, E, F), and three additional virtual sensors, which fall at the locations of other sensors (points B, C, E). The last three sensors, as discussed above, can be interpreted as spatial weighting or tapering. Finally, the virtual aperture of the orthogonal signals isAuthorized licensed use limited to: NORTHWESTERN POLYTECHNIC UNIVERSITY. Downloaded on January 7, 2010 at 03:53 from IEEE Xplore. Restrictions apply.BEKKERMAN AND TABRIKIAN: TARGET DETECTION AND LOCALIZATION USING MIMO RADARS AND SONARS3877created by nine sensors (six of them are at different locations), compared to three sensors of the coherent signals model. B. Spatial Coverage Extension In conventional target detection and localization systems, several directional beams are usually transmitted in order to scan a given region of interest (ROI). Each directional beam is generated using coherent transmitted signals. The time-on-target (TOT) for each transmitted beam is equal to the total interval assigned for covering the ROI, divided by the number of beams required to cover the given ROI. When the transmitted signals are orthogonal, the beams become omnidirectional, causing reduction in the beam gain. On the other hand, transmission of omnidirectional beams extends the spatial coverage of each beam; therefore, it allows increase of the TOT interval for each beam. In fact, the TOT interval for orthogonal signals is equal to the interval which is assigned to scan the ROI. Therefore, the beam gain loss can be compendirectional beams, one sated by increased TOT. Instead of omnidirectional beam can be transmitted with times higher TOT interval. Hence, in the sequel, when comparing the spatially orthogonal and coherent signals, this TOT compensation will be considered. Furthermore, the spatial transmitted power density with orthogonal signals is constant at each direction, while in the coherent signal case, the spatial transmitted power density is nonuniform and depends on the beam shape and the overlap between the beams. In omnidirectional signal transmission, the echo signal should be processed for a larger ROI, and therefore statistically a larger number of targets. If the targets in the given ROI are disjoint in range or Doppler, then for each range-Doppler bin, a single target case should be considered. In practice, the probability of multiple targets in a given range-Doppler bin is low, although it is still higher than the case of directional signal transmission in which only the targets within the narrow beam are excited. C. Beam Pattern Improvement The transmit–receive pattern can be written as (21) where is the target DOA and is the digital beam direction. Equation (21) can be rewritten in terms of the steering vector (see Appendix B) (22) It is worthwhile to notice that the right term in the numerator of (22) represents the beam pattern in the receive mode and is independent of the transmitted signal coherence matrix. The other terms in (22) represent the beam pattern in the transmit mode. For coherent transmit signals steered to the array broadside with , the transmit–receive pattern is given by (23)Fig. 5. Beam pattern for orthogonal and coherent signals with a ULA of 10 elements and with half a wavelength spacing.M=and for orthonormal signals creased TOT by factorconsidering the in-(24) in (23) introduces the attenuation in the The term transmit gain due to beam shape loss. This attenuation does not exist for the orthogonal signals case. The beam patterns for coherent and orthogonal transmitted signals are shown in Fig. 5. elements The array is uniform linear array (ULA) with and with half a wavelength spacing, where the transmit beam , i.e., the phases of of coherent signals is directed to . Fig. 5 shows that the orthogonal transmitted and lower signal model provides narrower beam width by sidelobe levels, compared to the coherent signal model. Note that the sidelobe level in the coherent signal case is at about 13 dB below the mainlobe level, which reﬂects the contribution of the digital beamforming in the receive mode only. In the orthogonal signal case, the sidelobe level is at about 26 dB below the mainlobe level reﬂecting the contribution of the digital beamforming in the transmission in addition to the receive mode. These phenomena can also be interpreted as the contribution of the virtual aperture extension and virtual tapering, as mentioned above. In Fig. 6, the beam patterns for coherent and orthogonal signals are shown, where the target is located at , , , respectively. It can be observed that for coherent transmitted signals, the gain of the transmit–receive pattern is attenuated when the target is not located in the center of the transmit beam, because of the beam shape loss. However, for orthogonal transmitted signals, the gain remains constant for all target directions . D. Increase of the Limit on the Number of Targets In conventional localization methods, an array of elements targets with unknown allows detection of up to complex amplitudes. In the spatially coded signal model, as discussed above, the number of different virtual sensors is theAuthorized licensed use limited to: NORTHWESTERN POLYTECHNIC UNIVERSITY. Downloaded on January 7, 2010 at 03:53 from IEEE Xplore. Restrictions apply.3878IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 54, NO. 10, OCTOBER 2006After optimization with respect to , the ML estimator for given byis(27) (28) where is a projection . For the matrix into the subspace spanned by columns of can be rewritten as scenario of a single target, (29) is a projection matrix into the subspace spanned by where . By substitution of from (12), from columns of (16), and using the identity (see [12]), with some trace properties, the numerator of (29) can be expressed as (30) is the sufﬁcient statistic matrix, deﬁned in (10). In where Appendix B, it is shown thatFig. 6. Beam pattern for orthogonal and coherent signals with a ULA of = 10 elements and with half a wavelength spacing. The beam in the transmit mode is directed to = 0 and the target is located at = 0 , 5 , 10 .M(31) Hence, (27) for a single target becomes (32)number of different combinations of transmit–receive delays according to (3). Hence, the limit on the maximum . number of targets is now constrained by This upper bound on the number of targets can be smaller depending on the array geometry. Thus, orthogonal signal transmission signiﬁcantly increases the upper limit on the number of targets that can be localized. V. TARGET DETECTION AND DOA ESTIMATION In order to demonstrate the advantages of the proposed conﬁguration in target detection and localization, the ML and CRB on DOA estimation and GLRT for target detection are derived.and the ML estimator ofcan be written as (33)B. Detection The hypotheses for a single target detection for the sufﬁcient statistic model in (15), can be stated as(34) Thus, the GLRT is given byA. Maximum Likelihood Estimation The model in (15) can be rewritten in matrix form(35) where and are the probability denand sity functions of the sufﬁcient statistic under hypotheses , respectively. Hence, the GLRT for the model in (6) can be written as (36) The threshold is set according to the desired false alarm rate. for It is interesting to ﬁnd the asymptotic statistics of(25) and . where Hence, the ML estimator for target localization for the model in (25) is given by(26)Authorized licensed use limited to: NORTHWESTERN POLYTECHNIC UNIVERSITY. Downloaded on January 7, 2010 at 03:53 from IEEE Xplore. Restrictions apply.BEKKERMAN AND TABRIKIAN: TARGET DETECTION AND LOCALIZATION USING MIMO RADARS AND SONARS3879both threshold setting and for performance evaluation. Unfor, the target parameters cannot tunately, under hypothesis be estimated, and thus derivation of the asymptotic properties of the above GLRT is a difﬁcult problem [13]. However, it can be veriﬁed that for orthogonal signals, there is no coupling between and . Therefore, the fact that is unknown does not affect the estimation performance of , which is required for decision between the two hypotheses. Accordingly, for the orunder thogonal signals case, the asymptotic statistics of the two hypotheses are given by [10] (37) where and are central and noncentral chi-squared distributions with two degrees of freedom, respectively, and is the noncentrality parameter, which is equal to (38) Using the Neyman–Pearson criterion, the probability of false alarm is const and the probability of detection is , where and are right-tail probability functions [10]. C. Cramér–Rao Bound The Fisher information matrix (FIM) [9] for estimating the vector from the sufﬁcient statistic , according to the model in (7), can be partitioned as (39) and the CRB for the DOA estimate can be expressed as CRB where , , and rived in Appendix C (40) CRB are matrices, whose elements are de. From (46), we where equality is satisﬁed only for conclude that the DOA estimation performance with orthogonal transmitted signals is superior to the performance obtained with , the TOT can coherent transmitted signals. Note that for , and thus the bound can be compensated by a factor of further be decreased. VI. SIMULATION RESULTS In Appendix C, it is shown that if the array origin is chosen to be at the array centroid, i.e., , the CRB for DOA estimation of a single target is given by (44) as shown at the . bottom of the page, where SNR For the case of two elements with (i.e., the transmit beam is steered to the array broadside), the CRB for In this section, we demonstrate via simulations the detection and localization performance for the case of spatially coded signals. A ULA with half a wavelength spacing is considered. In Fig. 7, the CRB for DOA estimation root-mean-square as a function of its error (RMSE) of a single target and orthonormal direction is plotted for coherent CRB (46)Fig. 7. CRB on DOA for = 10, L = 1, SNR= 0 dB. The target is located at = 0 and the beam is steered at for coherent signals.MDOA estimation can be expressed as a function of inserting CRB SNRby(45), , and is the where distance between the elements. The optimal value of , which minimizes the CRB, can be obtained by differentiating (45) with respect to and then equating to zero. The minimal CRB for is obtained for even without TOT compensation, which represents the case of orthonormal transmitted signals. Furthermore, it can be shown that(41) (42) (43)CRB SNR(44)Authorized licensed use limited to: NORTHWESTERN POLYTECHNIC UNIVERSITY. Downloaded on January 7, 2010 at 03:53 from IEEE Xplore. Restrictions apply.3880IEEE TRANSACTIONS ON SIGNAL PROCESSING, VOL. 54, NO. 10, OCTOBER 2006Fig. 8. CRB on DOA estimation, L = 2, M = 2, SNR= 0 dB; = 0 , = 5 , 10 , 15 . Two targets can be resolved by two elements if incoherent signals are transmitted.Fig. 10. ROC for the coherent and orthogonal transmitted signals; M = 10, L = 1, = 0 , and SNR= 10 dB. For the coherent transmitted signals, the beam is directed at = 0 .Fig. 9. CRB and ML performance for DOA estimation of the ﬁrst target for coherent and orthogonal transmitted signals; M = 10 elements, L = 2 targets: = 0 , 2 [0; 2] , SNR= 0 dB.signals . The array includes elements and dB. Note that for orthogonal transmitted signals, the SNR CRB is constant with respect to , while in the coherent transmitted signals case, the CRB increases with due to the beam shape loss. Fig. 8 presents the CRB for localization of a target located at in the presence of an additional target at , , with an array of elements. Obviously, the condition which is relevant to the case of coherent signals is not satisﬁed. However, Fig. 8 shows that the CRB is ﬁnite for . This phenomenon is directly related to the contribution of the virtual sensors in the noncoherent signals case, as mentioned in Section IV. In this case, the number ofvirtual sensors at the different locations is , and thus localization of two targets is possible with an array of two elements. As expected, the CRB goes to inﬁnity for . In Fig. 9, the angular resolution with coherent and orthogonal transmitted signals is examined. The scenario includes elements and targets, where the ﬁrst target is located at . The CRB and performance of the ML estimator of the ﬁrst target angle as a function of the position of the second for both coherent and orthogonal transmitted signals target are depicted. The TOT compensation is taken into consideration for the orthogonal transmitted signals. It can be noticed that the performance of the spatially coded signal model with orthogonal transmitted signals is superior to the conﬁguration with coherent transmitted signals.2 In Fig. 10, receiver operating characteristic (ROC) curves are presented using simulation results for both coherent and orthogonal transmitted signals. The theoretical asymptotic detection performance for the case of orthogonal signals is also calculated using (37) and depicted in this ﬁgure. The scenario inelements and target, which is located cludes at , and for the coherent transmitted signals the beam is . The TOT compensation is taken into considdirected at eration for the orthogonal transmitted signals. It can be observed that the detection performance with orthogonal signals is higher than with coherent signals. In addition, the asymptotic performance obtained theoretically coincides with the simulations for the case of orthogonal signals. VII. DISCUSSION AND CONCLUSIONS In this paper, a new approach for space-time signal transmission in radar and sonar systems was presented. Spatially orthogonal signal transmission enables digital beamforming and array2At low angular separation, the ML estimation error RMSE is lower than the CRB because of the limited search region for estimation of and .Authorized licensed use limited to: NORTHWESTERN POLYTECHNIC UNIVERSITY. Downloaded on January 7, 2010 at 03:53 from IEEE Xplore. Restrictions apply.。

胶州湾海水淤泥密度指数影响因素研究

DOI ：10.19965/ki.iwt.2022-1161第 43 卷第 11 期2023年 11 月Vol.43 No.11Nov.，2023工业水处理Industrial Water Treatment 胶州湾海水淤泥密度指数影响因素研究陈景光1，冯念林1，于崇涛2，杨鑫2，周利3（1.世帆工程设计有限公司，山东青岛 266034； 2.青岛海水淡化有限公司，山东青岛 266034；3.青岛理工大学环境与市政工程学院，山东青岛 266034）[ 摘要 ] 超滤+反渗透组合的双膜法在海水淡化中已经得到广泛应用，然而反渗透膜的污染问题严重影响了正常生产。

淤泥密度指数（SDI ）已经成为评价反渗透进水水质的主要指标。

通过对超滤出水水质进行检测分析并对反渗透膜主要污染物进行分析，探讨目前胶州湾海水淡化过程中反渗透进水SDI 的影响因素及各因素影响程度，结果表明：微生物作为自然水体中广泛存在的生物，其随季节和温度变化造成的总量变化是影响SDI 的最主要因素；水体中的硫酸盐以及Si 、Al 、Fe 胶体是造成膜污染的主要无机污染物；此外，浊度与水样的SDI 没有相关关系，也不能在反渗透过程中对水质情况进行有效衡量。

最后，针对影响SDI 的微生物因素的研究以及膜表面无机胶体的清洗提出未来发展方向。

［关键词］反渗透；膜污染；海水淡化；淤泥密度指数（SDI ）；微生物；胶体［中图分类号］ P747 ［文献标识码］B ［文章编号］ 1005-829X （2023）11-0208-07Influencing factors of seawater silt density index in Jiaozhou BayCHEN Jingguang 1，FENG Nianlin 1，YU Chongtao 2，YANG Xin 2，ZHOU Li 3（1.Shifan Engineering Design Co.， L td.， Q ingdao 266034，China ；2.Qingdao Seawater Desalination Co.， L td.， Q ingdao 266034，China ；3.School of Environmental and Municipal Engineering ，Qingdao University of Technology ，Qingdao 266034，China ）Abstract ：The combination of ultrafiltration and reverse osmosis has been widely used in seawater desalination ，but the pollution of reverse osmosis membrane has seriously affected the normal production. Silt density index （SDI ） has become the main index for evaluating reverse osmosis inlet water quality. By detecting and analyzing the water qual‑ity of ultrafiltration effluent and analyzing the main pollutants in the reverse osmosis membrane ，this study explored the influencing factors and the degree of influence of various factors on the SDI of reverse osmosis inflow seawater desalination process in Jiaozhou Bay. The results showed that ，as the widely existing organisms in natural water ，the total amount change of microorganisms with the season and temperature was the most important factor affecting SDI. Sulfate and Si ，Al ，Fe colloids in water were the main inorganic pollutants that caused membrane pollution. In addi‑tion ，turbidity had no correlation with SDI of water samples ，and could not be used to effectively measure water qual‑ity in reverse osmosis process. Finally ，future development directions were proposed for the study of microbial fac‑tors that affect SDI and the cleaning of inorganic colloids on membrane surfaces.Key words ：reverse osmosis ；membrane fouling ；seawater desalination ；silt density index （SDI ）；microorganisms ；colloids海水淡化即利用海水脱盐生产淡水，是实现水资源利用的开源增量技术，可以增加淡水总量，且不受时空和气候影响，保障沿海居民用水和工业用水等稳定供应。

空间句法的量化指标

空间句法的量化指标空间句法（Spatial Syntax）的量化指标主要用于描述和分析城市空间结构的特征和模式。

以下是一些常用的空间句法量化指标：1. 道路密度（Road Density）：指城市道路网络的密集程度，可通过计算单位面积内的道路长度来衡量。

2. 连通度（Connectivity）：指城市道路网络中的连通性程度，常用指标包括连接数（Connectivity Degree）、网络密度（Network Density）和平均最短路径长度（Average Shortest Path Length）等。

3. 可达性（Accessibility）：指城市不同区域的可达性程度，包括正向可达性（从某一区域到其他区域的便利程度）和反向可达性（从其他区域到某一区域的便利程度）。

4. 中心度（Centrality）：指城市空间中的中心程度，常用指标包括中心性（Centrality Degree）和中心距离（Centrality Distance）等。

5. 自由度（Integration）：指城市空间中的自由流动程度，包括自由度指数（Integration Index）和自由度梯度（Integration Gradient）等。

6. 分形维度（Fractal Dimension）：指城市空间形态的复杂程度，可通过计算城市边界、街道网络或建筑物形态的分形维度来衡量。

7. 线性维度（Linear Dimension）：指城市线性要素（如道路、河流等）的长度量化指标，包括总长度、平均长度和线性密度等。

8. 笑脸指数（Smile Index）：用于评估城市公共空间的多样性和流动性，从而衡量城市的活力和吸引力。

需要注意的是，以上仅列举了一些常见的空间句法量化指标，具体的指标选择和计算方法可能会因研究目的和所采用的空间分析工具而有所差异。

On the usage and measurement of landscape connectivity

OIKOS 90:7–19.Copenhagen 2000Minireviews provides an opportunity to summarize existing knowledge of selected ecological areas,with special emphasis on current topics where rapid and signiﬁcant advances are occurring.Reviews should be concise and not too wide-ranging.All key references should be cited.A summary is required.MINI-REVIEWOn the usage and measurement of landscape connectivityLutz Tischendorf and Lenore FahrigTischendorf,L.and Fahrig,L.2000.On the usage and measurement of landscape connectivity.–Oikos 90:7–19.This paper examines the usage and measurement of ‘‘landscape connectivity’’in 33recent studies.Connectivity is deﬁned as the degree to which a landscape facilitates or impedes movement of organisms among resource patches.However,connectivity is actually used in a variety of ways in the literature.This has led to confusion and lack of clarity related to (1)function vs structure,(2)patch isolation vs landscape connectivity and,(3)corridors vs connectivity.We suggest the term connectivity should be reserved for its original purpose.We highlight nine studies;these include modeling studies that actually measured connectivity in accordance with the deﬁni-tion,and empirical studies that measured key components of connectivity.We found that measurements of connectivity provide results that can be interpreted as recom-mending habitat fragmentation to enhance landscape connectivity.We discuss rea-sons for this misleading conclusion,and suggest a new way of quantifying connectivity,which avoids this problem.We also recommend a method for reducing sampling intensity in landscape-scale empirical studies of connectivity.L .Tischendorf and L .Fahrig ,Ottawa -Carleton Inst .of Biology ,Carleton Uni 6.,Ottawa ,ON ,Canada K 1S 5B 6(present address of LT :Busestrasse 76,D -28213Bremen ,Germany [tischendorf@cla 6is -bremen .de ]).What is landscape connectivity?The effects of spatial structure (patchiness)on popu-lation dynamics were ﬁrst examined in patch-based population models beginning in the early 1970s (e.g.,Levins 1969,Reddingius and den Boer 1970,Levin 1974,1976,Roff 1974).Further modeling studies showed that assumptions about movement among habitat patches greatly inﬂuence the predictions of such models (e.g.,Lefkovitch and Fahrig 1985,Fahrig 1988,1990,Fahrig and Paloheimo 1988,Henein and Merriam 1990,Adler and Nuernberger 1994,Lindenmayer and Lacy 1995,Lindenmayer and Possingham 1996,Frank and Wissel 1998,Henein et al.1998).Movement among habitat patches is,how-ever,not simply a function of an organism itself,but also depends on the landscape through which it must move.To emphasize the interaction between species’attributes and landscape structure in determining movement of organisms among habitat patches,Mer-riam (1984)introduced the concept of ‘‘landscape connectivity’’.OIKOS 90:1(2000)7Accepted 31January 2000Copyright ©OIKOS 2000ISSN 0030-1299Printed in Ireland –all rights reservedTaylor et al.(1993)deﬁned landscape connectivity as ‘‘the degree to which the landscape facilitates or im-pedes movement among resource patches’’.Similarly, With et al.(1997)deﬁned landscape connectivity as ‘‘the functional relationship among habitat patches, owing to the spatial contagion of habitat and the movement responses of organisms to landscape struc-ture’’.These deﬁnitions accentuate the dependence of movement on landscape structure,which suggests that connectivity is species-and landscape-speciﬁc.One must therefore describe landscape structure from a species’point of view(Wiens and Milne1989).This starts with deﬁning the species’habitat.The next step is to determine the scale at which the species responds to landscape structure,through itsﬁne-scale(grain)and large-scale(extent)movement(Wiens1997).This deter-mines the scale of habitat pattern as perceived by the organism.Finally,one must determine how the species responds to the different elements of a landscape.This comprises the species movement pattern and mortality risk on landscape elements(patches)as well as reactions at boundaries.Note that all of these behavioral facets contribute toward facilitating or impeding movement among resource patches.In summary,landscape connectivity encapsulates the combined effects of(1)landscape structure and(2)the species’use,ability to move and risk of mortality in the various landscape elements,on the movement rate among habitat patches in the landscape. Objective and approachWe reviewed the literature covered by the Agriculture, Biology&Environmental Sciences Edition of the Cur-rent Contents database(CC1998),from May1993to November1998.We searched article titles and key words for the term connecti6ity in combination with landscape or patch or habitat.The search resulted in49 papers.However,17of these papers did not use con-nectivity at all.We omitted these from the review,and included one other paper(Doak et al.1992)leaving33 papers,which are assembled in descending chronologi-cal and alphabetical order in Table1,and classiﬁed in Fig.1.Our objective was to examine the current usage and measurement of landscape connectivity.We start with a critical discussion of the diverse usage of connectivity, followed by a description of modeling and empirical studies that actually attempted to quantify connectivity or key components of it.We then discuss crucial model-ing assumptions and reveal the deceptive paradox of patch-based connectivity measurements,and its poten-tial for misleading conclusions.We end by suggesting ways to streamline and focus research on landscape connectivity.Current usage of connectivityStructure or function?The literature review revealed that the term connectiv-ity is sometimes used as a functional concept and other times in a structural way.Structural connectivity is equated with habitat contiguity and is measured by analyzing landscape structure,independent of any at-tributes of the organism(s)of interest(Collinge and Forman1998).The functional concept of connectivity explicitly con-siders the behavioral responses of an organism to the various landscape elements(patches and boundaries). Consequently,functional connectivity covers situations where organisms venture into non-habitat(matrix), where they may(1)face higher mortality risks(e.g., Lidicker1975,Gaines and McGlenaghan1980,Krohne and Burgin1987,Henein and Merriam1990,Schippers et al.1996,Charrier et al.1997,Poole1997,Sakai and Noon1997),(2)express different movement patterns (e.g.,Baars1979,Rijnsdorp1980,Wallin and Ekbom 1988,Wegner and Merriam1990,Hansson1991,John-son et al.1992a,Andreassen et al.1996b,Matter1996, Charrier et al.1997,Collins and Barrett1997),and(3) cross boundaries(e.g.,Mader1984,Wiens et al.1985, Bakowski and Kozakiewicz1988,Merriam et al.1989, Duelli et al.1990,Mader et al.1990,Frampton et al. 1995,Mauremooto et al.1995,Charrier et al.1997, Sakai and Noon1997).Depending on the movement attributes of the organ-ism,structural and functional connectivity can be syn-onymous.This occurs when the organism’s movement is conﬁned to its preferred habitat,i.e.,individuals do not cross the habitat/matrix boundary,and the organ-ism moves freely within the preferred habitat(e.g., Bascompte and Sole´1996).This is the assumption behind most percolation-based connectivity measures (Gardner et al.1987,Gardner and O’Neill1991,Green 1994).The fact that structural connectivity is relatively easy to measure could lead to the conclusion that connectiv-ity is a generalized feature of a landscape.This would be erroneous.In fact,the same landscape will have different connectivities for different organisms.Struc-turally connected habitat patches still may not be func-tionally connected and even non-contiguous habitat patches may be functionally connected,depending on the species(With1997).For example,if the only two habitat patches in a landscape are structurally con-nected by an inappropriate corridor for the species in question(too narrow or too long),structural connectiv-ity would exist without successful movement(functional response)from one patch to the other.Likewise,non-contiguous habitat patches may functionally be con-nected if the species can cross the non-habitat area (matrix)successfully and move between habitat8OIKOS90:1(2000)OIKOS 90:1(2000)9T a b l e 1.C h r o n o l o g i c a l a n d a l p h a b e t i c a l a s s e m b l a g e o f t h e 33r e v i e w e d c o n n e c t i v i t y s t u d i e s .S t u d y t y p e a n d d u r a t i o nS t u d yN o .M e a s u r e m e n t /u s a g e o f c o n n e c t i v i t yC o m m e n t s /s t u d y t a r g e t S p a t i a l s c a l e A n d r e a s s e n e t l a n d s c a p e e f f e c t s o n m o v e m e n t f r e q u e n c i e s e x p e r i m e n t ,r a d i o -t r a c k i n g ,13w k ,3p r e s e n c e /a b s e n c e o f c o r r i d o r s1g e n e r a t i o n s a l .1998e f f e c t o f s t r u c t u r a l p a t c h i s o l a t i o n o n s u m m e d i n ﬂu e n c e o f s i z e a n d s p a t i a l p a t c h o b s e r v a t i o n a l ,2y rA u l t a n d 2a r r a n g e m e n t o f n e i g hb o r i n g p a tc h e sc o m m u n i t y s t r u c t u r e a nd p o p u l a t i o n J o h n s o n 1998d e n s i t y e x p e r i m e n t ,l i v e -t r a p p i n g ,13w k ,3e f f e c t o n d i s p e r s a l d i s t a n c e s a n d s p a t i a l p r e s e n c e /a b s e n c e o f c o r r i d o r s l a n d s c a p eB j o r n s t a d e t a l .3a g g r e g a t i o n g e n e r a t i o n 1998l a n d s c a p e ,10e x p e r i m e n tp e r c e n t o f e q u a l l y s p a c e d s t r a i g h t l i n e s 4C o l l i n g e a n d e f f e c t o f s t r u c t u r a l c o n n e c t i v i t y m ×10mm e a s u r e m e n t o n i n s e c t d e n s i t y ,r i c h n e s s ,F o r m a n 1998c o v e r i n g h a b i t a t w i t h i n a l a nd s c a pe a n d c o m m u n i t y s t r u c t u r e d i s c u s s i o nr a t i o n a l e o n t h e q u a n t i ﬁc a t i o n o f l a n d s c a p es t r u c t u r a l p a t c h i s o l a t i o n D a v i d s o n 19985l a n d s c a p e f r a g m e n t a t i o n l a n d s c a p em o v e m e n t o f w a t e r b i r d s a m o n g h a b i t a t d i s c u s s i o nH a i g e t a l .1998r a t i o n a l e o n t h e i m p o r t a n c e o f f u n c t i o n a l 6p a t c h e s c o n n e c t i v i t y f o r w a t e r b i r d c o n s e r v a t i o n s p a t i a l l y e x p l i c i t s i m u l a t i o n m o d e l ,25y rH e n e i n e t a l .p r e s e n c e ,a b s e n c e a n d q u a l i t y o f f e n c e r o w s 7e f f e c t s o n p o p u l a t i o n s u r v i v a l l a n d s c a p ei n s i m u l a t e d l a n d s c a p e 1998e f f e c t o f s t r u c t u r a l a n d f u n c t i o n a l 8l a n d s c a p ee x p e r i m e n t ,o n e s e a s o nP e t i t a n d B u r e l d i s t a n c e s (e u c l i d i a n ,a l o n g h e d g e r o w s a n d ,w e i g h t e d b y m o v e m e n t i n t e n s i t y a n d c o n n e c t i v i t y o n s p e c i e s ’l o c a l a b u n d a n c e 1998bm o r t a l i t y i n d i f f e r e n t h a b i t a t t y p e s )b e t w e e n s a m p l e s i t e s 9l a n d s c a p ee x p e r i m e n t ,o n e s e a s o n‘‘f u n c t i o n a l d i s t a n c e ’’(w e i g h t e d b y P e t i t a n d B u r e l e f f e c t o f f u n c t i o n a l c o n n e c t i v i t y o n 1998am o v e m e n t i n t e n s i t y a n d m o r t a l i t y i n s p e c i e s ’l o c a l a b u n d a n c e d i f f e r e n t h a b i t a t t y p e s )l a n d s c a p em a n i p u l a t i v e m a r k -r e c a p t u r e e x p e r i m e n t ,e f f e c t s o f l a n d s c a p e s t r u c t u r e o n m o v e m e n t a b i l i t y o f d a m s e l ﬂi e s t h r o u g h 10P i t h e r a n d d i f f e r e n t h a b i t a t t y p e s m o v e m e n t f r e q u e n c i e s o n e s e a s o n T a y l o r 1998e f f e c t o f s t r u c t u r a l c o n n e c t i v i t y o n l a n d s c a p eG I S b a s e d p o p u l a t i o n d y n a m i c s m o d e l n e a r e s t n e i g h b o r d i s t a n c e (c o m b i n e d w i t h 11R o o t 1998d i s pe r s a lf r e q u e n c y d i s t r i b u t i o n )b e t w e e n (m e t a )p o p u l a t i o n s i z e (R A M A S )h a b i t a t p a t c h e s 12p a t c hv e g e t a t i o n s u r v e ya m o u n t o f f o r e s t h ab i t a t a r o u n d p a tc h e s e f f e c t o f s t r u c t u r a l p a t c h i s o l a t i o n o n G r a s h o f b o kd a m 1997w i t h i n t h r e e z o n e s u p t o 1000m z o o c h o r o u s a n d a n e m o c h o r o u s p l a n t s p e c i e s l a n d s c a p es t a t i c o p t i m i z a t i o n a n d s i m u l a t i o n m o d e lp o p u l a t i o n a b u n d a n c e r e l a t i o n s h i p 13H o f a n d c o n n e c t i v i t y m e a s u r e a s s u m e d t o b e b e t w e e n a d j a c e n t c e l l s i n a g r i d m o d e l s p a t i a l l i m i t a t i o n f a c t o r R a p h a e l 1997c o n n e c t i v i t y c o m p o n e n t s :a )d e g r e e o f r a t i o n a l e o n t h e e f f e c t o f c o n n e c t i v i t y o n l a n d s c a p ec o n c e p t u a l m ode l ,d i s c u s s i o n14M e t z g e r a n d h a b i t a t p e r c o l a t i o n ,b )c o r r i d o r a n d b i o d i v e r s i t y D e ´c a m p s 1997s t e p p i n g s t o n e n e t w o r k s ,c )m a t r i x p e r m e a b i l i t y e f f e c t o n l o c a l b i r d c o m m u n i t y (s p e c i e s S c h m i e g e l o w e t p a t c he x p e r i m e n t ,1y rc o r r id o r s (r i p a r i a n b u f fe r s t r i p s )b e t w e e n 15f o r e s t f r ag m e n t s a b u n d a n c e s )a l .1997e f f e c t s o f s i m u l a t e d l a n d s c a p e ch a n g e s o n p r o xi m i t y i n d e x -s u m m a r i z e d (p a t c h p a t c hq u a n t i ﬁc a t i o n o f s p a t i a l p a t t e r n i n G I SS p e t i c h e t a l .161997m a p ss t r u c t u r a l p a t c h i s o l a t i o n (p r o x i m i t y )a r e a /d i s t a n c e t o f o c a l p a t c h )r e l a t i o n s h i p f o r a l l p a t c h e s l o c a t e d w i t h i n r e c t a n g u l a r b u f f e r z o n e a r o u n d f o c a l p a t c h e s t i m a t i o n o f e f f e c t s o n s p e c i a l i s t l a n d s c a p ec o n c e p t u a l m ode li n t r i n s i c (j u x t a p o s i t i o n o f s i m i l a r h a b i t a t )T i e b o u t a n d 17c o l o n i z i n g a b i l i t y a nde x t r i n s i c (c o r r i d o r )c o n n e c t i v i t y A n d e r s o n 1997W i t h e t a l .ef f e c t o f l a n d s c a p e s p a t i a l s t r u c t u r e o n a v e r ag e d i s t a n c e b e t w e e n t w o s i t e s o f a l a n d s c a p e r a n d o m w a l k s i m u l a t i o n m o d e l o n 18n e u t r a l (r a n d o m a n d f r a c t a l )l a n d s c a p e 1997g r i d b e l o n g i n g t o t h e s a m e (p e r c o l a t i o n )p e r c o l a t i o n t h r e s h o l d a n d p o p u l a t i o n s ’s p a t i a l d i s t r i b u t i o nm a p sc l u s t e r ,p o p u l a t i o n s ’s p a t i a ld i s t r i b u t i o n10OIKOS 90:1(2000)T a b l e 1.(C o n t i n u e d )S p a t i a l s c a l e S t u d y t y p e a n d d u r a t i o nS t u d yN o .M e a s u r e m e n t /u s a g e o f c o n n e c t i v i t yC o m m e n t s /s t u d y t a r g e t p a t c h (c o r r i d o r )e x p e r i m e n t ,3m oA n d r e a s s e n e t a l .19s t r u c t u r a l d i s c o n t i n u i t i e s i n c o r r i d o r s e f f e c t s o n m o v e m e n t r a t e s 1996b e f f e c t s o n m o v e m e n t r a t e s p a t c h (c o r r i d o r )w i d t h o f c o r r i d o r s e x p e r i m e n t ,3m oA n d r e a s s e n e t a l .201996a l a n d s c a p e p h y s i c a l c o n n e c t i o n b e t w e e n p a t c h e s m e t a p o p u l a t i o n m o d e lH e s s 1996e f f e c t o n r e c o l o n i z a t i o n a n d e x t i n c t i o n 21(c o r r i d o r s )r a t e e f f e c t o f f u n c t i o n a l p a t c h i s o l a t i o n o n p a t c h h a b i t a t i s o l a t i o n b a s e d o n d i s p e r s a l o b s e r v a t i o n a l ,2y rH j e r m a n n a n d 22s p e c i e s p a t c h o c c u p a n c y d i s t a n c e d i s t r i b u t i o n s a n d n e g a t i v e I m s 1996e x p o n e n t i a l d i s p e r s a l f u n c t i o n e f f e c t o n p o p u l a t i o n m e a n s a n d 23p a t c h o p t i m i z a t i o n m o d e lb e t w e e n -p a tc h m o v e m e n t p r o b a b i l i t y H o f a nd F l a t he r v a r i a n c e s 1996d e p e n d e n t o n a )s p e c i e s d i s p e r s a l c a p a b i l i t y ,b )h a r s h n e s s of i n t e r -p a t c h e n v i r o n m e n t ,c )i n t e r -p a t c h d i s t a n c e e x p e r i m e n t ,1y rL e c o m t e a n d s i m u l a t e d p r e s e n c e /a b s e n c e o f 24e f f e c t s o n i n t e r -p a t c h d i s p e r s a l p a t c h C l o b e r t 1996c o r r i d o r s i n a n e x p e r i m e n t a l l a n d s c a p e p a t c h o b s e r v a t i o n a l ,1y rt o t a l l e n g t h o f h e d g e s i n a 0.5-k m P a i l l a t a n d B u t e t e f f e c t o f s t r u c t u r a l p a t c h i s o l a t i o n o n 25r a d i u s a r o u n d a s a m p l i n g p l o t (p a t c h )s p e c i e s a b u n d a n c e a n d ﬂu c t u a t i o n s 1996e f f e c t o f r e a l l a n d s c a p e s t r u c t u r e o n l a n d s c a p e G I S -r e l a t e d r a n d o m w a l k s i m u l a t i o n m o d e l26S c h i p p e r s e t a l .m o v e m e n t p r o b a b i l i t i e s b e t w e e n a l l f u n c t i o n a l c o n n e c t i v i t y p a i r s o f h a b i t a t p a t c h e s 1996G I S -r e l a t e d r a n d o m w a l k s i m u l a t i o n m o d e ld i s pe r s a l s u c c e s s r a t e ,f r a c t i o n o f S c h u m a k e r 1996e f f e c t o f l a n d s c a p e s t r u c t u r e o n 27l a n d s c a p e f u n c t i o n a l c o n n e c t i v i t y i n d i v i d u a l s t h a t l o c a t e d n e w t e r r i t o r i e s p r e s e n c e /a b s e n c e o f c o r r i d o r s i n e f f e c t o n m e t a p o p u l a t i o n p e r s i s t e n c e S w a r t a n d L a w e s l a n d s c a p e 28m o d e lm e t a p o p u l a t i o n m o d e l 1996G I S -r e l a t e d r a n d o m w a l k a n d p o p u l a t i o n d y n a m i c s 29D e m e r s e t a l .l a n d s c a p e e f f e c t o f r e a l l a n d s c a p e s t r u c t u r e o n d i s p e r s a l (c o l o n i z a t i o n )s u c c e s s s i m u l a t i o n m o d e l f u n c t i o n a l c o n n e c t i v i t y 1995l a n d s c a p e o b s e r v a t i o n a l ,3y rp r e s e n c e /a b s e n c e o f c o r r i d o r s a n d /o r 30e f f e c t o n m o v e m e n t r a t e s a n d d i s t a n c e s A r n o l d e t a l .1993s t e p p i n g s t o n e s s t a t i c o p t i m i z a t i o n m o d e ll a n d s c a p e o p t i m i z a t i o n o f h a b i t a t p l a c e m e n t p r o b a b i l i s t i c ,d i s t a n c e d e p e n d e n t H o f a n d J o y c e 31i s o l a t i o n o f a c e l l i n a g r i d m o d e l 1993d e g r e e t o w h i c h t h e l a n d s c a p e l a n d s c a p e c o n c e p t u a l d i s c u s s i o nT a y l o r e t a l .199332d e ﬁn i t i o n o f f u n c t i o n a l c o n n e c t i v i t y f a c i l i t a t e s o r i m p e d e s m o v e m e n t a m o n g r e s o u r c e p a t c h e s e f f e c t o f s c a l e o f c l u s t e r i n g o n s e a r c h t i m e -n u m b e r o f m o v e m e n t l a n d s c a p e r a n d o m w a l k s i m u l a t i o n m o d e l o n h i e r a r c h i c a l ,n e u t r a l 33D o a k e t a l .1992l a n d s c a p e m a p ss t e p s r e q u i r e d t o ﬁn d a n e w h a b i t a t f u n c t i o n a l c o n n e c t i v i t y p a t c hpatches.Research is needed to determine what,if any, simple measures of landscape structure can be used as measures of landscape connectivity.Patch isolation or landscape connectivity?Patch isolation is determined by the rate of immigration into the patch;the lower the immigration rate,the more isolated is the patch.Immigration rate depends on(1) the amount of occupied habitat surrounding the focal patch,(2)the number of emigrants leaving the sur-rounding habitat,(3)the nature of the intervening matrix,(4)the movement and perceptual abilities of the organism,and(5)the mortality risk of dispersers (Wiens et al.1993).Since(1)and(3)are landscape structural features and(4)and(5)are the organisms’responses to landscape structure,patch isolation de-pends on‘‘the degree to which the landscape facilitates or impedes movement...’’(Taylor et al.1993).Patch isolation is therefore imbedded within the concept of landscape connectivity.In fact,landscape connectivity is essentially equivalent to the inverse of the average degree of patch isolation over the landscape;a land-scape including mostly patches with a high degree of isolation will be less connected than vice versa.Five of the33studies we reviewed equated patch isolation with connectivity(Hjermann and Ims1996, Paillat and Butet1996,Grashofbokdam1997,Spetich et al.1997,Ault and Johnson1998).Even though patch isolation is clearly part of landscape connectivity (above),none of these studies estimated immigration rates into patches.Rather,they related a species’abun-dance or presence/absence in a patch to structural attributes of the surrounding landscape,such as dis-tance to the nearest occupied patch,or amount of habitat in a circle surrounding the patch.Such studies may reveal the relative importance of local patch vs surrounding landscape effects.However,they do not directly contribute to determining landscape connectiv-ity,because they do not actually determine rates of movement among patches.Corridors or connectivity?Corridors are narrow,continuous strips of habitat that structurally connect two otherwise non-contiguous habitat patches.The corridor concept(e.g.,Forman 1983,Bennett1990,Merriam1991,Saunders and Hobbs1991,Lindenmayer and Nix1993,Merriam and Saunders1993,Bonner1994,Dawson1994,Rosenberg et al.1997,Tischendorf1997a)originated from the generalized assumption that organisms do not venture into non-habitat.Under this assumption,addition of any habitat to a landscape increases the ability of organisms to move.Corridors in a landscape may therefore be a component of its connectivity if they promote movement among habitat patches,but they do not determine its connectivity.The degree to which corridors contribute to landscape connectivity depends on the nature of the corridors,the nature of the matrix and the response of the organism to both(Rosenberg et al.1997,Beier and Noss1998).Six of the reviewed studies equated the term connec-tivity with the presence/absence of corridors(Hess 1996,Lecomte and Clobert1996,Swart and Lawes 1996,Schmiegelow et al.1997,Andreassen et al.1998, Bjornstad et al.1998),and two studies associated con-nectivity with corridor width(Andreassen et al.1996a) or corridor continuity(Andreassen et al.1996b).The studies investigated(1)what features of a corridorFig.1.Classiﬁcation of the33reviewed studiesaccording to study type(a),year of publication(b),andusage of the termconnectivity(c).OIKOS90:1(2000)11determine its use by the organism,(2)space-use of organisms as a function of corridor presence/absence, and(3)population or community responses to corridors, e.g.,species richness,diversity or abundance.None of the studies explicitly recognized that corridors are only a component of the concept of landscape connectivity;they actually equated the connecting function of corridors with connectivity.Measurements of connectivityIn this section we review studies that quantiﬁed connec-tivity or key components of it.Recall that connectivity is deﬁned as the degree to which the landscape facilitates or impedes movement among resource patches.Only four of the studies(Doak et al.1992,Demers et al.1995, Schippers et al.1996,Schumaker1996)measured move-ments among resource patches over the entire landscape and actually quantiﬁed connectivity in accordance with its deﬁnition.All of these were modeling studies and were based on simulated movements across heterogeneous landscapes.We also reviewﬁve other studies which we think made an important contribution toward the con-cept of landscape connectivity(as explained below),even though they did not measure movement among resource patches directly(Arnold et al.1993,With et al.1997,Petit and Burel1998a,b,Pither and Taylor1998). Modeling studiesDispersal successDispersal success is usually deﬁned as the proportion of individuals that successfully immigrate into a new habitat patch during the course of a simulation run.Three of the modeling studies quantiﬁed connectivity using dispersal success.Schippers et al.(1996)(no.26in Table1)simulated the badger’s(Meles meles)response(movement proba-bility and mortality risk)to landscape structure using a classiﬁed GIS grid map and empirical expertise.Move-ment probabilities between cells were derived by compar-ing the quality(for badger use)of adjacent cells.Higher quality cells attracted moving individuals.Mortality rates were higher in low-quality cells.The number of simulated movement steps corresponded to an estimated actual time of badger movement within a four-year period.The authors produced inter-patch transition probabilities and movement frequency maps(visits per grid cell),based on dispersal success.Schumaker(1996)(no.27in Table1)analyzed the potential of indices of landscape structure to predict dispersal success.He created landscape models in two ways:(1)sample landscapes were randomly drawn from a GIS data set to cover a range of different landscape conﬁgurations;(2)artiﬁcial landscape grids were created by randomly designating habitat cells.Cells of the grid represented territories.An individual-based correlated random walk model was used to simulate movements across the landscape.Individuals were released in a randomly selected50%of habitat territories,and were allowed to settle in any unoccupied territory,which then became unavailable to subsequent nd-scape boundaries reﬂected approaching individuals. Connectivity was measured as the mean fraction(over several runs)of individuals that successfully dispersed into new territories during the course of a simulation.The results revealed correlations between each of ten indices of landscape structure and dispersal success(connectiv-ity).Demers et al.(1995)(no.29in Table1)investigated the relationship between colonization success of edge-preferring organisms,and the amount and change of edge habitat,in real agricultural landscapes.A vector-based GIS data set containing fencerow and forest-edge cover-ages was used as a model landscape.Individuals were allowed to move only in suitable habitat after being dropped at random points across the landscape.Individ-uals could cross inhospitable habitat(matrix)up to a maximum distance,after any edge habitat in the land-scape was successfully colonized.Occupied habitat could not be colonized by subsequent dispersers.The authors measured connectivity as the‘‘total length and area of hedgerow and forest edge colonized by the offspring of each successful virtual organism’’.The results showed higher connectivity in landscapes with more and longer overall edge habitat.Search timeOne paper(Doak et al.1992)(no.33in Table1)used search time to quantify connectivity.Search time is the number of movement steps individuals require toﬁnd a new habitat patch.Doak et al.(1992)examined the effect of spatial scale on the success of dispersing individuals.An artiﬁcial landscape was modeled by a hierarchical grid of three layers(spatial scales).Clusters of habitat cells were created on different spatial scales.Virtual individuals were released in the habitat and followed a random walk until a new habitat patch(different from the origin)was ndscape boundaries acted as reﬂecting borders. For each individual the number of movement steps required toﬁnd a new habitat patch(search time)was recorded.The mean and standard deviation over all individual search times were calculated and related to the scale of rge-scale clustering(few large patches)induced longer search times than small-scale clustering(more smaller patches)(see also Ruckelshaus et al.1997).Population spatial distributionWith et al.(1997)(no.18in Table1)investigated the effects of landscape spatial structure on(1)the probabil-12OIKOS90:1(2000)。

Advances in

Advances in Geosciences,4,17–22,2005 SRef-ID:1680-7359/adgeo/2005-4-17 European Geosciences Union©2005Author(s).This work is licensed under a Creative CommonsLicense.Advances in GeosciencesIncorporating level set methods in Geographical Information Systems(GIS)for land-surface process modelingD.PullarGeography Planning and Architecture,The University of Queensland,Brisbane QLD4072,Australia Received:1August2004–Revised:1November2004–Accepted:15November2004–Published:9August2005nd-surface processes include a broad class of models that operate at a landscape scale.Current modelling approaches tend to be specialised towards one type of pro-cess,yet it is the interaction of processes that is increasing seen as important to obtain a more integrated approach to land management.This paper presents a technique and a tool that may be applied generically to landscape processes. The technique tracks moving interfaces across landscapes for processes such as waterﬂow,biochemical diffusion,and plant dispersal.Its theoretical development applies a La-grangian approach to motion over a Eulerian grid space by tracking quantities across a landscape as an evolving front. An algorithm for this technique,called level set method,is implemented in a geographical information system(GIS).It ﬁts with aﬁeld data model in GIS and is implemented as operators in map algebra.The paper describes an implemen-tation of the level set methods in a map algebra program-ming language,called MapScript,and gives example pro-gram scripts for applications in ecology and hydrology.1IntroductionOver the past decade there has been an explosion in the ap-plication of models to solve environmental issues.Many of these models are speciﬁc to one physical process and of-ten require expert knowledge to use.Increasingly generic modeling frameworks are being sought to provide analyti-cal tools to examine and resolve complex environmental and natural resource problems.These systems consider a vari-ety of land condition characteristics,interactions and driv-ing physical processes.Variables accounted for include cli-mate,topography,soils,geology,land cover,vegetation and hydro-geography(Moore et al.,1993).Physical interactions include processes for climatology,hydrology,topographic landsurface/sub-surfaceﬂuxes and biological/ecological sys-Correspondence to:D.Pullar(d.pullar@.au)tems(Sklar and Costanza,1991).Progress has been made in linking model-speciﬁc systems with tools used by environ-mental managers,for instance geographical information sys-tems(GIS).While this approach,commonly referred to as loose coupling,provides a practical solution it still does not improve the scientiﬁc foundation of these models nor their integration with other models and related systems,such as decision support systems(Argent,2003).The alternative ap-proach is for tightly coupled systems which build functional-ity into a system or interface to domain libraries from which a user may build custom solutions using a macro language or program scripts.The approach supports integrated models through interface speciﬁcations which articulate the funda-mental assumptions and simpliﬁcations within these models. The problem is that there are no environmental modelling systems which are widely used by engineers and scientists that offer this level of interoperability,and the more com-monly used GIS systems do not currently support space and time representations and operations suitable for modelling environmental processes(Burrough,1998)(Sui and Magio, 1999).Providing a generic environmental modeling framework for practical environmental issues is challenging.It does not exist now despite an overwhelming demand because there are deep technical challenges to build integrated modeling frameworks in a scientiﬁcally rigorous manner.It is this chal-lenge this research addresses.1.1Background for ApproachThe paper describes a generic environmental modeling lan-guage integrated with a Geographical Information System (GIS)which supports spatial-temporal operators to model physical interactions occurring in two ways.The trivial case where interactions are isolated to a location,and the more common and complex case where interactions propa-gate spatially across landscape surfaces.The programming language has a strong theoretical and algorithmic basis.The-oretically,it assumes a Eulerian representation of state space,Fig.1.Shows a)a propagating interface parameterised by differ-ential equations,b)interface fronts have variable intensity and may expand or contract based onﬁeld gradients and driving process. but propagates quantities across landscapes using Lagrangian equations of motion.In physics,a Lagrangian view focuses on how a quantity(water volume or particle)moves through space,whereas an Eulerian view focuses on a localﬁxed area of space and accounts for quantities moving through it.The beneﬁt of this approach is that an Eulerian perspective is em-inently suited to representing the variation of environmen-tal phenomena across space,but it is difﬁcult to conceptu-alise solutions for the equations of motion and has compu-tational drawbacks(Press et al.,1992).On the other hand, the Lagrangian view is often not favoured because it requires a global solution that makes it difﬁcult to account for local variations,but has the advantage of solving equations of mo-tion in an intuitive and numerically direct way.The research will address this dilemma by adopting a novel approach from the image processing discipline that uses a Lagrangian ap-proach over an Eulerian grid.The approach,called level set methods,provides an efﬁcient algorithm for modeling a natural advancing front in a host of settings(Sethian,1999). The reason the method works well over other approaches is that the advancing front is described by equations of motion (Lagrangian view),but computationally the front propagates over a vectorﬁeld(Eulerian view).Hence,we have a very generic way to describe the motion of quantities,but can ex-plicitly solve their advancing properties locally as propagat-ing zones.The research work will adapt this technique for modeling the motion of environmental variables across time and space.Speciﬁcally,it will add new data models and op-erators to a geographical information system(GIS)for envi-ronmental modeling.This is considered to be a signiﬁcant research imperative in spatial information science and tech-nology(Goodchild,2001).The main focus of this paper is to evaluate if the level set method(Sethian,1999)can:–provide a theoretically and empirically supportable methodology for modeling a range of integral landscape processes,–provide an algorithmic solution that is not sensitive to process timing,is computationally stable and efﬁcient as compared to conventional explicit solutions to diffu-sive processes models,–be developed as part of a generic modelling language in GIS to express integrated models for natural resource and environmental problems?The outline for the paper is as follow.The next section will describe the theory for spatial-temporal processing us-ing level sets.Section3describes how this is implemented in a map algebra programming language.Two application examples are given–an ecological and a hydrological ex-ample–to demonstrate the use of operators for computing reactive-diffusive interactions in landscapes.Section4sum-marises the contribution of this research.2Theory2.1IntroductionLevel set methods(Sethian,1999)have been applied in a large collection of applications including,physics,chemistry,ﬂuid dynamics,combustion,material science,fabrication of microelectronics,and computer vision.Level set methods compute an advancing interface using an Eulerian grid and the Lagrangian equations of motion.They are similar to cost distance modeling used in GIS(Burroughs and McDonnell, 1998)in that they compute the spread of a variable across space,but the motion is based upon partial differential equa-tions related to the physical process.The advancement of the interface is computed through time along a spatial gradient, and it may expand or contract in its extent.See Fig.1.2.2TheoryThe advantage of the level set method is that it models mo-tion along a state-space gradient.Level set methods start with the equation of motion,i.e.an advancing front with velocity F is characterised by an arrival surface T(x,y).Note that F is a velocityﬁeld in a spatial sense.If F was constant this would result in an expanding series of circular fronts,but for different values in a velocityﬁeld the front will have a more contorted appearance as shown in Fig.1b.The motion of thisinterface is always normal to the interface boundary,and its progress is regulated by several factors:F=f(L,G,I)(1)where L=local properties that determine the shape of advanc-ing front,G=global properties related to governing forces for its motion,I=independent properties that regulate and inﬂu-ence the motion.If the advancing front is modeled strictly in terms of the movement of entity particles,then a straightfor-ward velocity equation describes its motion:|∇T|F=1given T0=0(2) where the arrival function T(x,y)is a travel cost surface,and T0is the initial position of the interface.Instead we use level sets to describe the interface as a complex function.The level set functionφis an evolving front consistent with the under-lying viscosity solution deﬁned by partial differential equa-tions.This is expressed by the equation:φt+F|∇φ|=0givenφ(x,y,t=0)(3)whereφt is a complex interface function over time period 0..n,i.e.φ(x,y,t)=t0..tn,∇φis the spatial and temporal derivatives for viscosity equations.The Eulerian view over a spatial domain imposes a discretisation of space,i.e.the raster grid,which records changes in value z.Hence,the level set function becomesφ(x,y,z,t)to describe an evolv-ing surface over time.Further details are given in Sethian (1999)along with efﬁcient algorithms.The next section de-scribes the integration of the level set methods with GIS.3Map algebra modelling3.1Map algebraSpatial models are written in a map algebra programming language.Map algebra is a function-oriented language that operates on four implicit spatial data types:point,neighbour-hood,zonal and whole landscape surfaces.Surfaces are typ-ically represented as a discrete raster where a point is a cell, a neighbourhood is a kernel centred on a cell,and zones are groups of mon examples of raster data include ter-rain models,categorical land cover maps,and scalar temper-ature surfaces.Map algebra is used to program many types of landscape models ranging from land suitability models to mineral exploration in the geosciences(Burrough and Mc-Donnell,1998;Bonham-Carter,1994).The syntax for map algebra follows a mathematical style with statements expressed as equations.These equations use operators to manipulate spatial data types for point and neighbourhoods.Expressions that manipulate a raster sur-face may use a global operation or alternatively iterate over the cells in a raster.For instance the GRID map algebra (Gao et al.,1993)deﬁnes an iteration construct,called do-cell,to apply equations on a cell-by-cell basis.This is triv-ially performed on columns and rows in a clockwork manner. However,for environmental phenomena there aresituations Fig.2.Spatial processing orders for raster.where the order of computations has a special signiﬁcance. For instance,processes that involve spreading or transport acting along environmental gradients within the landscape. Therefore special control needs to be exercised on the order of execution.Burrough(1998)describes two extra control mechanisms for diffusion and directed topology.Figure2 shows the three principle types of processing orders,and they are:–row scan order governed by the clockwork lattice struc-ture,–spread order governed by the spreading or scattering ofa material from a more concentrated region,–ﬂow order governed by advection which is the transport of a material due to velocity.Our implementation of map algebra,called MapScript (Pullar,2001),includes a special iteration construct that sup-ports these processing orders.MapScript is a lightweight lan-guage for processing raster-based GIS data using map alge-bra.The language parser and engine are built as a software component to interoperate with the IDRISI GIS(Eastman, 1997).MapScript is built in C++with a class hierarchy based upon a value type.Variants for value types include numeri-cal,boolean,template,cells,or a grid.MapScript supports combinations of these data types within equations with basic arithmetic and relational comparison operators.Algebra op-erations on templates typically result in an aggregate value assigned to a cell(Pullar,2001);this is similar to the con-volution integral in image algebras(Ritter et al.,1990).The language supports iteration to execute a block of statements in three ways:a)docell construct to process raster in a row scan order,b)dospread construct to process raster in a spreadwhile(time<100)dospreadpop=pop+(diffuse(kernel*pop))pop=pop+(r*pop*dt*(1-(pop/K)) enddoendwhere the diffusive constant is stored in thekernel:Fig.3.Map algebra script and convolution kernel for population dispersion.The variable pop is a raster,r,K and D are constants, dt is the model time step,and the kernel is a3×3template.It is assumed a time step is deﬁned and the script is run in a simulation. Theﬁrst line contained in the nested cell processing construct(i.e. dospread)is the diffusive term and the second line is the population growth term.order,c)doﬂow to process raster byﬂow order.Examples are given in subsequent sections.Process models will also involve a timing loop which may be handled as a general while(<condition>)..end construct in MapScript where the condition expression includes a system time variable.This time variable is used in a speciﬁc fashion along with a system time step by certain operators,namely diffuse()andﬂuxﬂow() described in the next section,to model diffusion and advec-tion as a time evolving front.The evolving front represents quantities such as vegetation growth or surface runoff.3.2Ecological exampleThis section presents an ecological example based upon plant dispersal in a landscape.The population of a species follows a controlled growth rate and at the same time spreads across landscapes.The theory of the rate of spread of an organism is given in Tilman and Kareiva(1997).The area occupied by a species grows log-linear with time.This may be modelled by coupling a spatial diffusion term with an exponential pop-ulation growth term;the combination produces the familiar reaction-diffusion model.A simple growth population model is used where the reac-tion term considers one population controlled by births and mortalities is:dN dt =r·N1−NK(4)where N is the size of the population,r is the rate of change of population given in terms of the difference between birth and mortality rates,and K is the carrying capacity.Further dis-cussion of population models can be found in Jrgensen and Bendoricchio(2001).The diffusive term spreads a quantity through space at a speciﬁed rate:dudt=Dd2udx2(5) where u is the quantity which in our case is population size, and D is the diffusive coefﬁcient.The model is operated as a coupled computation.Over a discretized space,or raster,the diffusive term is estimated using a numerical scheme(Press et al.,1992).The distance over which diffusion takes place in time step dt is minimally constrained by the raster resolution.For a stable computa-tional process the following condition must be satisﬁed:2Ddtdx2≤1(6) This basically states that to account for the diffusive pro-cess,the term2D·dx is less than the velocity of the advancing front.This would not be difﬁcult to compute if D is constant, but is problematic if D is variable with respect to landscape conditions.This problem may be overcome by progressing along a diffusive front over the discrete raster based upon distance rather than being constrained by the cell resolution.The pro-cessing and diffusive operator is implemented in a map al-gebra programming language.The code fragment in Fig.3 shows a map algebra script for a single time step for the cou-pled reactive-diffusion model for population growth.The operator of interest in the script shown in Fig.3is the diffuse operator.It is assumed that the script is run with a given time step.The operator uses a system time step which is computed to balance the effect of process errors with efﬁ-cient computation.With knowledge of the time step the it-erative construct applies an appropriate distance propagation such that the condition in Eq.(3)is not violated.The level set algorithm(Sethian,1999)is used to do this in a stable and accurate way.As a diffusive front propagates through the raster,a cost distance kernel assigns the proper time to each raster cell.The time assigned to the cell corresponds to the minimal cost it takes to reach that cell.Hence cell pro-cessing is controlled by propagating the kernel outward at a speed adaptive to the local context rather than meeting an arbitrary global constraint.3.3Hydrological exampleThis section presents a hydrological example based upon sur-face dispersal of excess rainfall across the terrain.The move-ment of water is described by the continuity equation:∂h∂t=e t−∇·q t(7) where h is the water depth(m),e t is the rainfall excess(m/s), q is the discharge(m/hr)at time t.Discharge is assumed to have steady uniformﬂow conditions,and is determined by Manning’s equation:q t=v t h t=1nh5/3ts1/2(8)putation of current cell(x+ x,t,t+ ).where q t is theﬂow velocity(m/s),h t is water depth,and s is the surface slope(m/m).An explicit method of calcula-tion is used to compute velocity and depth over raster cells, and equations are solved at each time step.A conservative form of aﬁnite difference method solves for q t in Eq.(5). To simplify discussions we describe quasi-one-dimensional equations for theﬂow problem.The actual numerical com-putations are normally performed on an Eulerian grid(Julien et al.,1995).Finite-element approximations are made to solve the above partial differential equations for the one-dimensional case ofﬂow along a strip of unit width.This leads to a cou-pled model with one term to maintain the continuity ofﬂow and another term to compute theﬂow.In addition,all calcu-lations must progress from an uphill cell to the down slope cell.This is implemented in map algebra by a iteration con-struct,called doﬂow,which processes a raster byﬂow order. Flow distance is measured in cell size x per unit length. One strip is processed during a time interval t(Fig.4).The conservative solution for the continuity term using aﬁrst or-der approximation for Eq.(5)is derived as:h x+ x,t+ t=h x+ x,t−q x+ x,t−q x,txt(9)where the inﬂow q x,t and outﬂow q x+x,t are calculated in the second term using Equation6as:q x,t=v x,t·h t(10) The calculations approximate discharge from previous time interval.Discharge is dynamically determined within the continuity equation by water depth.The rate of change in state variables for Equation6needs to satisfy a stability condition where v· t/ x≤1to maintain numerical stabil-ity.The physical interpretation of this is that aﬁnite volume of water wouldﬂow across and out of a cell within the time step t.Typically the cell resolution isﬁxed for the raster, and adjusting the time step requires restarting the simulation while(time<120)doflow(dem)fvel=1/n*pow(depth,m)*sqrt(grade)depth=depth+(depth*fluxflow(fvel)) enddoendFig.5.Map algebra script for excess rainfallﬂow computed over a 120minute event.The variables depth and grade are rasters,fvel is theﬂow velocity,n and m are constants in Manning’s equation.It is assumed a time step is deﬁned and the script is run in a simulation. Theﬁrst line in the nested cell processing(i.e.doﬂow)computes theﬂow velocity and the second line computes the change in depth from the previous value plus any net change(inﬂow–outﬂow)due to velocityﬂux across the cell.cycle.Flow velocities change dramatically over the course of a storm event,and it is problematic to set an appropriate time step which is efﬁcient and yields a stable result.The hydrological model has been implemented in a map algebra programming language Pullar(2003).To overcome the problem mentioned above we have added high level oper-ators to compute theﬂow as an advancing front over a land-scape.The time step advances this front adaptively across the landscape based upon theﬂow velocity.The level set algorithm(Sethian,1999)is used to do this in a stable and accurate way.The map algebra script is given in Fig.5.The important operator is theﬂuxﬂow operator.It computes the advancing front for waterﬂow across a DEM by hydrologi-cal principles,and computes the local drainageﬂux rate for each cell.Theﬂux rate is used to compute the net change in a cell in terms ofﬂow depth over an adaptive time step.4ConclusionsThe paper has described an approach to extend the function-ality of tightly coupled environmental models in GIS(Ar-gent,2004).A long standing criticism of GIS has been its in-ability to handle dynamic spatial models.Other researchers have also addressed this issue(Burrough,1998).The con-tribution of this paper is to describe how level set methods are:i)an appropriate scientiﬁc basis,and ii)able to perform stable time-space computations for modelling landscape pro-cesses.The level set method provides the following beneﬁts:–it more directly models motion of spatial phenomena and may handle both expanding and contracting inter-faces,–is based upon differential equations related to the spatial dynamics of physical processes.Despite the potential for using level set methods in GIS and land-surface process modeling,there are no commercial or research systems that use this mercial sys-tems such as GRID(Gao et al.,1993),and research systems such as PCRaster(Wesseling et al.,1996)offerﬂexible andpowerful map algebra programming languages.But opera-tions that involve reaction-diffusive processing are speciﬁc to one context,such as groundwaterﬂow.We believe the level set method offers a more generic approach that allows a user to programﬂow and diffusive landscape processes for a variety of application contexts.We have shown that it pro-vides an appropriate theoretical underpinning and may be ef-ﬁciently implemented in a GIS.We have demonstrated its application for two landscape processes–albeit relatively simple examples–but these may be extended to deal with more complex and dynamic circumstances.The validation for improved environmental modeling tools ultimately rests in their uptake and usage by scientists and engineers.The tool may be accessed from the web site .au/projects/mapscript/(version with enhancements available April2005)for use with IDRSIS GIS(Eastman,1997)and in the future with ArcGIS. It is hoped that a larger community of users will make use of the methodology and implementation for a variety of environmental modeling applications.Edited by:P.Krause,S.Kralisch,and W.Fl¨u gelReviewed by:anonymous refereesReferencesArgent,R.:An Overview of Model Integration for Environmental Applications,Environmental Modelling and Software,19,219–234,2004.Bonham-Carter,G.F.:Geographic Information Systems for Geo-scientists,Elsevier Science Inc.,New York,1994. Burrough,P.A.:Dynamic Modelling and Geocomputation,in: Geocomputation:A Primer,edited by:Longley,P.A.,et al., Wiley,England,165-191,1998.Burrough,P.A.and McDonnell,R.:Principles of Geographic In-formation Systems,Oxford University Press,New York,1998. Gao,P.,Zhan,C.,and Menon,S.:An Overview of Cell-Based Mod-eling with GIS,in:Environmental Modeling with GIS,edited by: Goodchild,M.F.,et al.,Oxford University Press,325–331,1993.Goodchild,M.:A Geographer Looks at Spatial Information Theory, in:COSIT–Spatial Information Theory,edited by:Goos,G., Hertmanis,J.,and van Leeuwen,J.,LNCS2205,1–13,2001.Jørgensen,S.and Bendoricchio,G.:Fundamentals of Ecological Modelling,Elsevier,New York,2001.Julien,P.Y.,Saghaﬁan,B.,and Ogden,F.:Raster-Based Hydro-logic Modelling of Spatially-Varied Surface Runoff,Water Re-sources Bulletin,31(3),523–536,1995.Moore,I.D.,Turner,A.,Wilson,J.,Jenson,S.,and Band,L.:GIS and Land-Surface-Subsurface Process Modeling,in:Environ-mental Modeling with GIS,edited by:Goodchild,M.F.,et al., Oxford University Press,New York,1993.Press,W.,Flannery,B.,Teukolsky,S.,and Vetterling,W.:Numeri-cal Recipes in C:The Art of Scientiﬁc Computing,2nd Ed.Cam-bridge University Press,Cambridge,1992.Pullar,D.:MapScript:A Map Algebra Programming Language Incorporating Neighborhood Analysis,GeoInformatica,5(2), 145–163,2001.Pullar,D.:Simulation Modelling Applied To Runoff Modelling Us-ing MapScript,Transactions in GIS,7(2),267–283,2003. Ritter,G.,Wilson,J.,and Davidson,J.:Image Algebra:An Overview,Computer Vision,Graphics,and Image Processing, 4,297–331,1990.Sethian,J.A.:Level Set Methods and Fast Marching Methods, Cambridge University Press,Cambridge,1999.Sklar,F.H.and Costanza,R.:The Development of Dynamic Spa-tial Models for Landscape Ecology:A Review and Progress,in: Quantitative Methods in Ecology,Springer-Verlag,New York, 239–288,1991.Sui,D.and R.Maggio:Integrating GIS with Hydrological Mod-eling:Practices,Problems,and Prospects,Computers,Environ-ment and Urban Systems,23(1),33–51,1999.Tilman,D.and P.Kareiva:Spatial Ecology:The Role of Space in Population Dynamics and Interspeciﬁc Interactions.Princeton University Press,Princeton,New Jersey,USA,1997. Wesseling C.G.,Karssenberg, D.,Burrough,P. A.,and van Deursen,W.P.:Integrating Dynamic Environmental Models in GIS:The Development of a Dynamic Modelling Language, Transactions in GIS,1(1),40–48,1996.。

朱章志运用扶正祛邪法论治糖尿病经验

ʌ临证验案ɔ朱章志运用扶正祛邪法论治糖尿病经验❋曾绘域1,朱章志2ә,周㊀海3,陈㊀珺3,张文婧3(1.深圳市中西医结合医院,广东深圳㊀518104;2.广州中医药大学第一附属医院,广州㊀510405;3.广州中医药大学,广州㊀510405)㊀㊀摘要:糖尿病属于中医学消渴病范畴,以往医家多认为其病机为阴虚燥热,治疗以滋阴清热为法㊂朱章志教授通过长期的临床观察与实践,立足于张仲景保胃气,扶阳气的理论,认为糖尿病的病机为正虚邪滞,即太阴虚损㊁阳气不足㊁收敛不及,寒㊁水㊁湿之邪阻滞阳气运行通道㊂治疗上不囿陈法,以扶正祛邪为大法,通过固护太阴㊁扶助阳气㊁收敛阳气,祛除寒水湿之邪,恢复阳气运行之通畅,使阳气功能复常㊁运行有序,为糖尿病的治疗提供临床新思路㊂㊀㊀关键词:扶正祛邪;糖尿病;朱章志㊀㊀中图分类号:R587.1㊀㊀文献标识码:A㊀㊀文章编号:1006-3250(2021)01-0149-03Discussion on ZHU Zhang-zhi's Experience in Treating Diabetes Mellitus by Using The Method of Reinforcing The Healthy Qi and Eliminating The Pathogenic FactorsZENG Hui-yu 1,ZHU Zhang-zhi 2ә,ZHOU Hai 3,CHEN Jun 3,ZHANG Wen-jing 3(1.Shenzhen Hospital of Integrated traditional Chinese and Western Medicine,Guangdong,Shenzhen 518104,China;2.The First Affiliated Hospital of Guangzhou University of Chinese Medicine,Guangzhou 510405,China;3.Guangzhou University of Chinese Medicine,Guangzhou 510405,China)㊀㊀Abstract :Diabetes mellitus belongs to the category of "xiao ke"in traditional Chinese medicine.Doctors used to think that its pathogenesis was Yin deficiency and dryness heat ,and the treatment was nourishing Yin and clearing heat.Through long-term clinical observation and practice ,and based on ZHANG Zhong-jing's theory of protecting stomach Qi and supporting Yng Q ,professor ZHU Zhang-zhi believes that the pathogenesis of diabetes is deficiency of healthy Qi and stagnation of pathogen.Because of the deficiency of greater Yin and Yang Qi ,and the lack of convergence ,the cold ,water and dampness block the operational channel of Yang Qi.The treatment of diabetes mellitus should be based on reinforcing the healthy Qi and eliminating the pathogenic factors.By strengthening Taiyin ,supporting Yang Qi ,astringent Yang Qi ,dispelling the evil of cold ,water and dampness ,we can restore the smooth operation of Yang Qi ,restore the function of Yang Qi to normal and operate orderly ,which provides a new clinical method for the treatment of diabetes mellitus.㊀㊀Key words :Reinforcing the healthy Qi and eliminating the pathogenic factors ;Diabetes mellitus ;ZHU Zhang-zhi❋基金项目:国家自然科学基金资助项目(81873190)-降糖三黄片在糖脂毒性所致胰岛β细胞损伤的自噬调控作用作者简介:曾绘域(1990-),女,广东云浮人,住院医师,硕士研究生,从事六经辨治内分泌疾病的临床与研究㊂ә通讯作者:朱章志(1963-),男,湖南衡阳人,主任医师,博士研究生导师,从事六经辨治内分泌疾病的临床与研究,Tel :************,E-mail :zhuangi@ ㊂㊀㊀随着人口老龄化和生活方式的改变,我国糖尿病的患病率呈上升趋势,2013年我国18岁以上人群糖尿病患病率为10.4%[1]㊂中医药在延缓糖尿病的进展及防治其并发症方面具有一定优势[2-4]㊂糖尿病属于中医学消渴病范畴,以往医家多认为其病机为阴虚燥热,治疗以滋阴清热为法,但疗效尚不能令人满意㊂朱章志教授通过长期的临床观察与实践,认为正虚邪滞乃糖尿病病机之核心,采用扶正祛邪法治之屡获奇效㊂1㊀正虚邪滞之糖尿病病机‘素问㊃经脉别论篇“曰: 饮入于胃,游溢精气,上输于脾,脾气散精水精四布,五经并行㊂食物入胃,经脾胃运化化生精气,然后输布全身㊂糖尿病患者常嗜食肥甘,起居无常,烦劳紧张,致太阴虚损,正气内虚,阳气戕伐,津液代谢异常,而生寒水湿之邪㊂寒㊁水㊁湿之邪气作为阴邪,又可阻滞阳气运行之通道㊂阳气运行通道不畅,不能敷布温煦四肢,可见手足逆冷;阳气运行受阻,又可出现郁而化热之象㊂因此朱章志认为,疗糖尿病的关键在于恢复阳气运行之通畅,根据糖尿病正虚邪滞的病机,治疗以扶正祛邪为法,顾护太阴㊁扶助阳气㊁收敛阳气,祛除寒水湿之邪,使阳气功能复常则行有序㊂2㊀运用扶正祛邪法治疗糖尿病2.1㊀扶正2.1.1㊀固护中气,扶助阳气㊀张仲景遣方用药常体现保胃气之思想[5],如桂枝汤中配伍生姜㊁大枣㊁炙甘草,发汗祛邪不忘顾护中气;又如白虎汤中加梗米㊁炙甘草以和中益胃,又可防止石膏㊁知母大寒伤中㊂有胃气则生,无胃气则死 ,故扶正之要以保胃气为先㊂朱章志认为,阳气在人体的生命活动中占主导9412021年1月第27卷第1期January 2021Vol.27.No.1㊀㊀㊀㊀㊀㊀中国中医基础医学杂志Journal of Basic Chinese Medicine地位㊂‘素问㊃生气通天论篇“曰: 阳气者若天与日,失其所则折寿而不彰是故阳因而上,卫外者也㊂ ‘黄帝内经“把阳气比作太阳,阳气运行失常可致短寿㊂阳气具有抵御外邪㊁护卫生命㊁维持机体生命活动的作用,津液的气化㊁血液的运行均需阳气的温煦与推动㊂因此,在人体的阴阳平衡中阳气起着主导作用㊂朱章志认为,正气虚衰㊁太阴虚损㊁阳气不足是糖尿病发生发展之根本原因,因此扶正首当固护中气,扶助阳气 ,故常以附子理中汤为底方,固护中宫㊂太阴脾土居中央,犹如足球比赛之中场,能联系前锋与后卫进可攻退可守,进可充养肺卫之气抵御外邪,退可顾护少阴以防寒邪内陷㊂‘四圣心源㊃卷二太阴湿土“提到: 湿者,太阴土气之所化也故胃家之燥,不敌脾家之湿,病则土燥者少而土湿者多也㊂[6] 阴脾土易挟寒湿,附子理中汤功善固护中气㊁温补脾阳而散寒湿,为治疗太阴阳虚寒湿之要方㊂方中附子辛温大热,补坎中真阳,又能散寒湿,荡去群阴;干姜去脏腑沉寒痼冷,温暖脾土,复兴火种;人参被誉为百草之王能大补元气,为扶正固本之极品;白术味苦性温,功善健脾燥湿,乃扶植太阴之要药;炙甘草善益气补中,调和药性,诸药合用以收培补中阳㊁散寒除湿之效㊂若其人神疲懒言,气虚较甚,在附子理中汤的基础上可重用红参㊁北芪以大补元气,健脾益气;若其人四肢不温㊁肢体困重㊁寒湿较重者,可加重附子㊁干姜之量,并加细辛㊁吴茱萸以散久寒;若其人口干口苦㊁舌苔黄腻㊁大便黏滞不爽兼夹湿热之象,可仿当归拈痛汤之意,加茵陈㊁当归㊁黄芩以利湿清热㊂2.1.2㊀收敛阳气,阳密乃固㊀朱章志认为, 阴可理解为阳气的收敛㊁收藏状态,糖尿病阴虚燥热之象乃阳气不足㊁收敛不及㊁升发太过所致[7]㊂‘素问㊃生气通天论篇“提到: 阳气者,烦劳则张㊂现代人起居无节,以妄为常,阳气因而不能潜藏,常常浮越于外容易出现假热之象,医者不察,妄投清热泻火之品,实乃雪上加霜㊂凡阴阳之要,阳密乃固 ,收敛阳气即是扶正,犹如太极之能收能放,收敛是为了聚集能量,阳气固密,正气才能强盛,方能更好的制敌㊂朱章志常用砂仁㊁肉桂㊁白芍㊁山萸肉㊁泽泻等药物收敛阳气㊂砂仁辛温,既能宣太阴之寒湿,又能纳气归肾,使阳气收敛于少阴,少火生气㊂‘本草经疏“提到: 缩砂蜜,辛能散,又能润辛以润肾,故使气下行气下则气得归元㊂[8] 肉桂引火归原,导浮越之阳气归于命门,益火消阴㊂若患者出现咽痛㊁牙龈肿痛㊁痤疮等阳气不敛㊁虚火上冲之象,常用砂仁㊁肉桂以收敛阳气,纳气归肾,引火归原㊂白芍味酸能敛,敛降甲木胆火,使相火归位㊂‘本草求真“曰: 气之盛者,必赖酸为之收,故白芍号为敛肝之液,收肝之气,而令气不妄行也㊂[9] 朱章志常使用白芍以补肝之体㊁助肝之用,收敛肝气,肝平则郁气自除,火热自消㊂山萸肉秘精气㊁敛阳气,使龙雷之火归于水中㊂朱章志常用山萸肉收敛正气,遇汗出多者,常重用以固涩敛汗㊂泽泻能泻能降,能入肾泻浊,开气化之源,泻浊以利扶正,又能降气而引火下行㊂朱章志常用泽泻打通西方潜藏之要塞[10],在温阳之品中加入泽泻,利于阳气潜藏,使孤阳有归㊂2.1.3㊀填补阴精,以滋化源㊀‘素问㊃金匮真言论篇“提到: 夫精者,身之本也㊂精是人体生命活动的物质基础,能化气生髓,濡养脏腑㊂人体之精禀受于父母,又由后天水谷之精不断充养,归藏于肾中㊂孤阴不生,独阳不长 ,无阳则阴无以生,无阴则阳无以化㊂肾乃水火之脏,阴精充足才能涵养坎中真火,使真阳固密于内,化生正气㊂朱章志常在秋冬之季嘱糖尿病患者进补阿胶等血肉有情之品填补肾精㊂肾主封藏,秋冬进补使肾精充养,以滋阳气化生之源㊂阿胶用黄酒烊化,既能祛除阿胶之腥,又能借黄酒通行之性解阿胶滋腻碍胃之弊,每日少量服用,以有形之精难以速生,填补肾精以缓补为要㊂除此之外,遣方用药时亦会注意顾护阴精,在使用温阳药的同时常常配伍山萸肉㊁白芍等养阴药,以防温燥伤阴之弊㊂2.2㊀祛邪2.2.1㊀外散寒水以运太阳㊀太阳为开 ,太阳乃三阳之表,巨阳也,其性开泄以应天,为祛邪之重要通道㊂在运气里,太阳在天为寒,在地为水,合而为太阳寒水㊂张仲景太阳病篇研究的是水循环过程,治太阳就是治水[11]㊂寒㊁水之邪闭郁在表,气血运行不畅,可见肌肤麻木不仁㊂邪气滞留太阳,阻碍阳气运行,当因势利导㊁开太阳之表以发汗,外散寒㊁水之邪㊂糖尿病患者正气亏虚为本,祛邪不能伤正,朱章志临床常用桂枝麻黄各半汤小发其汗,使玄府开张,邪有出路㊂桂枝麻黄各半汤乃发汗轻剂,为桂枝汤与麻黄汤相合而得,其中麻黄㊁桂枝㊁生姜㊁北杏发散宣肺以开皮毛,芍药㊁大枣㊁炙甘草酸甘化阴以益营,诸药相合,刚柔相济,祛邪而不伤正㊂邪去正安,阳气运行通畅,水液代谢复常则阳气自充,而无寒水之扰㊂若寒邪较重可用三拗汤,此为麻黄汤去桂枝而成,功善开宣肺气,疏散风寒,因去辛温之桂枝发汗力不及麻黄汤,祛邪而不伤正㊂2.2.2㊀下利水湿以健少阴㊀少阴乃水火交会之脏,元气之根,人身立命之本㊂‘医理真传“提到: 坎中真阳,一名龙雷火,发而为病,一名元阳外越,一名孤阳上浮,一名虚火上冲㊂此际之龙,乃初生之龙,不能飞腾而兴云布雨,惟潜于渊中,以水为家,以水为性,遂安其在下之位㊂水盛一分龙亦盛一分,水高一尺龙亦高一尺,是龙之因水盛而游 [12]㊂阴盛051中国中医基础医学杂志Journal of Basic Chinese Medicine㊀㊀㊀㊀㊀㊀2021年1月第27卷第1期January2021Vol.27.No.1则阳衰,水湿之邪泛滥,则龙雷之火因而飞越在外㊂叶天士深谙张仲景之理,提到通阳不在温,而在利小便 [10,13],通过利小便的方法,使水湿之邪从下而解,阳气运行通道无水湿之邪阻碍,则阳气无需温养而自通,水盛得除则真龙亦安其位㊂朱章志常用五苓散㊁真武汤下利水湿,以复阳气之通达,少阴之健运㊂五苓散具有通阳化气利水之效,治疗膀胱气化不利形成的蓄水证㊂方中猪苓㊁茯苓㊁泽泻导水湿之邪下行;白术健脾燥湿,杜绝生湿之源;桂枝助膀胱气化,通阳化气行水又通气于表,使全身在表之湿邪皆得解,五药合用,膀胱气化复常,水道通调使小便得利,水湿得出㊂真武汤为治疗少阴阳虚㊁水气泛滥之主方,方中附子振奋少阴阳气,使水有所主;白术㊁茯苓健脾制水;生姜助附子敷布阳气,宣散水气;芍药利小便,制附㊁姜之燥,五味相合共奏温阳利水之功㊂2.2.3㊀开郁逐寒以畅厥阴㊀肝为将军之官,肝气主动主升发,能统帅兵马,捍卫君主㊂厥阴肝经,体阴用阳,内寄相火,相火敷布阳气,祛阴除寒,是祛邪的先锋主力军㊂朱章志常用吴茱萸汤祛除厥阴肝经之寒邪,恢复肝经阳气之运行㊂方中吴茱萸辛苦而温,芳香而燥,‘本草汇言“曰: 开郁化滞,逐冷降气之药 [14],能温胃暖肝,降浊阴止呕逆,为治疗肝寒之要药㊂配以生姜温胃散寒,佐以人参㊁大枣健脾益气补虚,全方散寒与降逆并施,共奏暖肝温胃㊁降逆止呕之效㊂‘素问㊃至真要大论篇“说: 帝曰:厥阴何也?岐伯曰:两阴交尽也㊂物极必反,重阴必阳㊂厥阴为阴尽阳生之脏,足厥阴肝经与足少阳胆经互为表里,若出现肝气不疏㊁枢机不利㊁气郁化火,朱章志常用小柴胡汤和畅枢机,开郁以复气机调达㊂方中柴胡疏泄肝胆之气;黄芩清泄胆火,一疏一清,气郁通达,火郁得发;生姜㊁半夏和胃降逆;人参㊁大枣㊁炙甘草固护中宫,全方寒温并用㊁补泻兼施,以复厥阴疏泄之职,使气机和畅㊁阳气运行有序㊂3㊀典型病案患者杨某,女,65岁,2017年3月10日初诊:2型糖尿病病史6年余,症见疲乏,双下肢轻度浮肿,下肢冰凉,背部易汗出,口苦口干,偶有腰膝酸软,纳眠可,二便调,舌淡暗,苔黄腻,脉沉细㊂辅助检查示糖化血红蛋白10.8%,空腹血糖14.59mmol/L,总胆固醇6.38mmol/L,甘油三酯3.66mmol/L,低密度脂蛋白胆固醇4.34mmol/L㊂西医诊断2型糖尿病㊁高脂血症,治疗给予门冬胰岛素30(早餐前22u晚餐前20u)控制血糖,阿托伐他汀钙片(20mg, qn)调脂㊂中医诊断消渴病,少阴阳虚寒湿证㊂患者少阴阳气衰微不足以养神,固见疲乏;腰为肾之府,少阴阳虚则见腰膝酸软,阳虚寒盛则见下肢冰凉;背部正中乃督脉运行之所,阳气虚衰无以固摄则见背部汗出;少阴阳虚不能主水,寒水泛滥则见双下肢浮肿;水湿内停有郁而化热之象,则见口苦口干㊁舌苔黄腻㊁舌淡暗,脉沉细为少阴阳虚寒湿之征,治以温阳散寒㊁利水除湿为法㊂方以扶正祛邪方合当归拈痛汤加减:炮附片10g(先煎1h),红参10g (另炖),干姜15g,白术30g,炙甘草15g,桂枝12 g,麻黄8g,生姜30g,猪苓10g,泽泻30g,茯苓30 g,白芍30g,酒萸肉45g,当归15g,茵陈10g,5剂水煎服,2d1剂,水煎至250ml,饭后分2次服用,次日复煎㊂方中以附子理中汤为底方温补中焦,散寒除湿;加桂枝㊁麻黄使寒湿之邪从皮毛而解;加五苓散通阳化气,使湿邪从下而出;生姜散寒除湿;白芍㊁酒萸肉收敛阳气,以助正气祛邪;当归活血利水;茵陈清热利湿㊂2017年3月24日二诊:患者双下肢浮肿减轻,疲乏较前好转,无口干口苦,无腰膝酸软,仍觉下肢冰凉,背部仍有汗出,动则尤甚,大便每日二行,质偏烂,舌淡暗,苔白腻,脉细㊂患者大便质烂,乃邪有出路,导水湿之邪从大便而解㊂患者无口干口苦,舌苔由黄腻转为白腻,知湿郁化热之象已除,遂去茵陈㊂仍觉下肢冰凉乃内有久寒,加制吴茱萸12g以散沉寒痼冷;上方加酒萸肉至60g以加强收敛阳气㊁固摄敛汗之效,加黄芪60g以健脾益气敛汗;加砂仁6g(后下)㊁肉桂3g(焗服)以加强收敛阳气㊁扶助正气之效,7剂水煎服,服法同前㊂2017年4月7日三诊:患者背部汗出减少,下肢转温,余症皆除,大便每日二行质软,舌淡红,苔薄白,脉细较前有力,继续服二诊方药5剂㊂后给予附子理中丸(每次8粒,每日3次)服用1个月巩固疗效㊂2017年11月17日复诊:患者上述症状皆除,纳眠可,二便调㊂复查糖化血红蛋白6.8%,空腹血糖6.5mmol/L,总胆固醇5.14mmol/L,甘油三酯1.65 mmol/L,低密度脂蛋白胆固醇2.43mmol/L㊂4㊀结语以往医家多以滋阴清热为法治疗糖尿病,通过长期的临床实践,朱章志不囿陈法,根据糖尿病患者当前之病因病机特点,运用扶正祛邪法治疗糖尿病,通过顾护太阴㊁扶助阳气㊁收敛阳气,祛除寒水湿之邪气,恢复阳气运行之通畅,为糖尿病的治疗提供新思路㊂参考文献:[1]㊀WANG L GAO-P-ZHANG-M,et al.Prevalence and EthnicPattern of Diabetes and Prediabetes in China in2013[J].JAMA,2017,317(24):2515-2523.[2]㊀谭宏韬,刘树林,朱章志,等. 首辨阴阳,再辨六经论治惠州地区2型糖尿病的临床观察[J].中华中医药杂志,2018,33(9):4240-4244.(下转第181页)offspring of sleep-deprived mice[J].Psychoneuroendocrinology,2009,35(5):775-784.[9]㊀覃甘梅,覃骊兰.心肾不交型失眠动物模型研究进展[J].中华中医药杂志,2018,33(1):229-231.[10]㊀吕志平,刘承才.肝郁致瘀机理探讨[J].中医杂志,2000,41(6):367-368.[11]㊀游秋云,王平,田代志,等.老年肝郁失眠证候大鼠模型的建立和评价[J].中国实验方剂学杂志,2013,19(2):222-225. [12]㊀唐仕欢,杨洪军,黄璐琦. 以方测证方法应用的反思[J].中国中西医结合杂志,2007,27(3):259-262.[13]㊀卢岩,刘振华,于晓华,等.疏肝调神针法针刺对睡眠剥夺模型大鼠神经递质的影响[J].山东中医杂志,2017,36(4):322-325. [14]㊀YANG CR,SEAMANS JK,GORELOVA N.Developing aneuronal model for the pathophysiology of schizophrenia based onthe nature of electrophysiological actions of dopamine in theprefrontal cortex[J].Neuropsychopharmacology,1999,21(2):161-194.[15]㊀何林熹,诸毅晖,杨翠花,等.失眠肝郁化火证大鼠模型的建立及其评价[J].中华中医药杂志,2018,33(9):3890-3894. [16]㊀李越峰,徐富菊,张泽国,等.四逆散对大鼠睡眠时相影响的实验研究[J].中国临床药理学杂志,2014,30(10):936-938. [17]㊀张晓婷,刘文超,刘俊昌,等.电击法建立SD大鼠焦虑型心理应激-失眠模型的研究[J].现代中西医结合杂志,2018,27(30):3316-3319.[18]㊀钱伯初,史红,郑晓亮.新的失眠动物模型研究概述[J].中国新药杂志,2008,17(1):1-4.[19]㊀朱洁,申国明,汪远金,等.肝郁证失眠大鼠模型的建立与评价[J].中医杂志,2011,52(8):689-692.[20]㊀刘倩,李蜀平,廖磊,等.调和肝脾方治疗失眠的实验研究[J].北京中医药,2018,37(8):768-770.[21]㊀全世建,焦蒙蒙,黑赏艳,等.交泰丸对睡眠剥夺大鼠下丘脑Orexin A及γ-氨基丁酸的影响[J].广州中医药大学学报,2015,32(1):103-105.[22]㊀KOBAN M,SWINSON KL.Chronic REM-sleep deprivation ofrats elevates metabolic rate and increases UCP1gene expressionin brown adipose tissue[J].Am J Physiol Endocrinol Metab,2005,289(1):68-74.[23]㊀赵俊云,杨晓敏,胡秀华,等.失眠动物模型HPA轴和表观遗传修饰的变化及交泰丸的干预作用[J].中医药学报,2018,46(4):19-21.[24]㊀GORGULU Y,CALIYURT O.Rapid antidepressant effects ofsleep deprivation therapy correlates with serum BDNF changes inmajor depression[J].Brain Res Bull,2009,80(3):158-162.[25]㊀BENCA RM,PETERSON MJ.Insomnia and depression[J].Sleep Med,2008,9(1):S3-S9.[26]㊀郜红利,涂星,卢映,等.心肾不交所致失眠大鼠模型[J].中成药,2014,36(6):1138-1141.[27]㊀杨钰涵,孙雨,王珺,等.中医病证相符的大鼠心肾不交失眠模型的建立及其血清代谢组学研究[J].中国中药杂志,2020,45(2):383-390.[28]㊀石皓月,鲁艺,李钰昕,等.中药治疗对氯苯丙氨酸失眠模型大鼠影响的基础研究进展[J].中国医药导报,2018,15(11):33-36.[29]㊀全世建,何树茂,钱莉莉.交泰丸交通心肾治疗失眠作用机理研究[J].辽宁中医药大学学报,2011,13(8):12-14. [30]㊀GULEC M,OZKOL H,SELVI Y,et al.Oxidative stress inpatients with primary insomnia[J].Pro NeuropsychopharmacolBiol Psychiatry,2012,37(2):247-251.[31]㊀ZHANG H,CAO D,CUI W,et al.Molecular bases ofthioredoxin and thioredoxin reductase-mediated prooxidant actionsof(-)-epigallocatechin-3-gallate[J].Free Radic Biol Med,2010,49(12):2010-2018.[32]㊀谢光璟,刘源才,胡辉,等.基于Trx系统介导的抗氧化应激探讨天王补心方对失眠模型大鼠的干预作用[J].时珍国医国药,2019,30(4):805-808.[33]㊀黄攀攀,王平,李贵海,等.老年阴虚失眠动物模型的建立与评价[J].中华中医药学刊,2010,28(8):1719-1723.[34]㊀XIONG L,HUANG XJ,SONG PX.The experiment ofstudymodel of Deficiency of yin Insomnia by Yangyin anshenkoufuye[J].Chin J Pract Chin Mod Med,2005,18(18):1187-1188.[35]㊀韦祎,唐汉庆,李克明,等.脾阳虚证失眠大鼠模型的建立和附子理中汤的干预效应[J].中国实验方剂学杂志,2013,19(16):289-292.[36]㊀王志鹏.桂枝甘草龙骨牡蛎汤对阳虚证失眠大鼠脑内5-HT㊁NE含量的影响[D].南京:东南大学,2015.[37]㊀MURRAY NM,BUCHANAN GF,RICHERSON GB.InsomniaCaused by Serotonin Depletion is Due to Hypothermia[J].Sleep,2015,38(12):1985-1993.[38]㊀宋亚刚,李艳,崔琳琳,等.中医药病证结合动物模型的现代应用研究及思考[J].中草药,2019,50(16):3971-3978. [39]㊀李晓娟,白晓晖,陈家旭,等.中医动物模型研制方法及展望[J].中华中医药杂志,2014,29(7):2263-2266.[40]㊀刘臻,谢晨,赵娜,等.失眠动物模型的制作与评价[J].中医学报,2013,28(12):1846-1848.收稿日期:2020-05-18(上接第151页)[3]㊀司芹芹,牛晓红,杨海卿,等.温阳益气养阴活血方治疗2型糖尿病肾病的临床疗效[J].中华中医药学刊,2018,36(3):703-705.[4]㊀郭仪,石岩.中药复方治疗糖尿病大血管病变临床疗效及对血糖㊁血脂影响的系统评价[J].中华中医药学刊,2017,35(6):1369-1375.[5]㊀方春平,刘步平,朱章志.‘伤寒论“中胃气思想在病脉辨证中的运用[J].浙江中医药大学学报,2014,38(8):948-950.[6]㊀黄元御.四圣心源[M].北京:中国中医药出版社,2009:24.[7]㊀朱章志,林明欣,樊毓运.立足阳主阴从浅析糖尿病的中医治疗[J].江苏中医药,2011,43(4):7-8.[8]㊀缪希雍.本草经疏[M].北京:中国医药科技出版社,2011:56.[9]㊀黄宫绣.本草求真[M].北京:中国中医药出版社,2008:132.[10]㊀林明欣,裴倩,朱章志.浅析通阳不在温,而在利小便 [J].中医杂志,2011,52(19):1705-1706.[11]㊀刘力红.思考中医[M].桂林:广西师范大学出版社,2006:457.[12]㊀郑寿全.医理真传[M].北京:中国中医药出版社,2008:3.[13]㊀刘涛,张毅,李娟,等.结合‘伤寒论“探讨通阳不在温而在利小便 [J].中国中医药信息杂志,2017,24(9):106-107.[14]㊀倪朱谟.本草汇言[M].北京:中医古籍出版社,2010:87.收稿日期:2020-04-27。

英语导游试题及答案

英语导游试题及答案一、选择题（每题2分，共20分）1. The most widely spoken language in the world is ________.A. EnglishB. MandarinC. SpanishD. French答案：A2. Which of the following is NOT a tourist attraction in London?A. The British MuseumB. The LouvreC. The Tower of LondonD. Big Ben答案：B3. The term "check-in" refers to ________.A. registering at a hotelB. boarding an airplaneC. checking out of a hotelD. leaving an airplane答案：A4. The phrase "see the sights" is commonly used to mean________.A. to visit tourist attractionsB. to watch a movieC. to look at the sceneryD. to meet with someone答案：A5. What is the meaning of the acronym "UNESCO"?A. United Nations Educational, Scientific and Cultural OrganizationB. United Nations Economic and Social CouncilC. United Nations Environmental and Social CommitteeD. United Nations Emergency Service Corps答案：A6. The phrase "hop on" is often used to describe ________.A. getting off a vehicleB. boarding a vehicleC. speeding up a vehicleD. stopping a vehicle答案：B7. In the context of tourism, "package tour" refers to________.A. a tour that includes all expensesB. a tour that is organized by a travel agencyC. a tour that is self-guidedD. a tour that focuses on outdoor activities答案：A8. The term "cultural exchange" involves ________.A. exchanging money for different currenciesB. sharing and learning about different culturesC. exchanging gifts between countriesD. trading goods between different countries答案：B9. What does "off the beaten track" mean in the context of travel?A. a popular tourist destinationB. a place that is not commonly visited by touristsC. a place that is difficult to reachD. a place that is on the main road答案：B10. The phrase "to go sightseeing" is synonymous with________.A. going shoppingB. going to the cinemaC. visiting places of interestD. going to the beach答案：C二、填空题（每题1分，共10分）11. The ________ is the official language of the UnitedNations.答案：English12. The Eiffel Tower is located in ________.答案：Paris13. A tourist who is interested in art might visit the________ in Florence, Italy.答案：Uffizi Gallery14. The term "budget travel" refers to traveling with a________.答案：limited amount of money15. A traveler might use a ________ to find the best route to their destination.答案：map16. The Great Barrier Reef is a popular destination for________.答案：scuba diving17. A ________ is a person who guides tourists and provides information about places of interest.答案：tour guide18. The ________ is a famous landmark in New York City.答案：Statue of Liberty19. A traveler might need a ________ to enter a foreign country.答案：visa20. The Colosseum is an ancient amphitheater located in________.答案：Rome三、简答题（每题5分，共30分）21. What are the responsibilities of a tour guide?答案：A tour guide's responsibilities include providing information about tourist attractions, ensuring the safety of tourists, and facilitating a smooth and enjoyable travel experience.22. What are some common challenges faced by tourists when traveling abroad?答案：Some common challenges include language barriers, cultural differences, unfamiliarity with local customs, and potential issues with transportation or accommodations.23. Explain the concept of "eco-tourism."答案：Eco-tourism is a form of tourism that focuses on sustainability and minimal impact on the environment. It often involves visiting natural areas and learning about conservation efforts.24. What is the significance of learning a few phrases in the local language when traveling?答案：Learning a few phrases in the local language can help tourists communicate more effectively, show respect for the local culture, and enhance their overall travel experience.四、论述题（每题15分，共30分）25. Discuss the importance of cultural sensitivity when working as an English-speaking tour guide in a。

安徽宏村古村落非物质文化遗产数字化展示研究

安徽建筑大学学报第30卷104水平较低，对于智慧城市的规划较少，相对来说评价指标较少。

此外，智慧城建评价体系不仅应根据本地经济发展、建设水平的提高进行修改，更应根据宏观行业发展规划的变化进行适当调整，以满足不断变化的实际需求，如城市住建行业已经与城市发展、居民日常生活密不可分，建筑节能也已经成为安徽省必不可少的节能政策之一。

因此，合肥、省辖市及区县均涉及相关的智慧城市评价指标。

各地根据智慧城建评价指标体系对智慧城建的发展建设情况进行考核，以达到促进行业发展，提高行业竞争力，带动本地相关产业以及经济总体的发展，更好地为建设智慧城市，提升居民幸福感，促进城市和谐、健康发展服务。

参考文献：[1] 李德仁，姚远，邵振峰.智慧城市的概念、支撑技术及应用[J].工程研究-跨学科视野中的工程，2012，4（4）：313-323.[2] Gil-Garcia J R，Pardo T A，Nam T. What makes a city smart? Identifying core components and proposing an integrative and comprehensive conceptualization[J].Information Polity，2015，20（1）：61-87.[3] Khatoun R，Zeadally S. Smart cities[J].Communications of the ACM，2016，59（8）：46-57.[4] Meijer A，Bolívar M P R. Governing the smart city：a review of the literature on smart urban governance[J].International Review of Administrative Sciences，2016， 82（2）：392-408.[5] 马惠雯.大数据背景下智慧城市建设的创新路径[J].中小企业管理与科技（中旬刊），2021（11）：46-48. [6] 王朝南.智慧城市建设绩效评价指标体系构建及实证研究[D].湘潭：湘潭大学，2019.[7] 安永文.张掖市智慧城市建设研究[D].武汉：湖北工业大学，2018.[8] 中国智慧城市建设行业现状研究分析及市场前景预测报告（2016年）[R].北京：中国产业调研网，2016. [9] 陈铭，王乾晨，张晓海，等.“智慧城市”评价指标体系研究——以“智慧南京”建设为例[J].城市发展研究，2011，18（5）：84-89.[10] 陈桂龙.智慧城市2.0的“浦东模式”[J].中国建设信息，2015（13）：34-36.[11] 王思雪，郑磊.国内外智慧城市评价指标体系比较[J].电子政务，2013（1）：92-100.[12] 常文辉.智慧城市评价指标体系构建研究[D].开封：河南大学，2014.[13] 单志广.国家新型智慧城市评价数据分析报告[R].沈阳：国家信息中心，2017.[14] 丁国胜，宋彦.智慧城市与“智慧规划”——智慧城市视野下城乡规划展开研究的概念框架与关键领域探讨[J].城市发展研究，2013，21（8）：34-39.[15] 郭素娴.智慧城市评价指标体系的构建及应用[D].杭州：浙江工商大学，2013.[16] 安徽省统计局.安徽统计年鉴（2018汉英对照附光盘）[M].北京：中国统计出版社，2018.[28] Yu Z，Hu J. Microstructure and characteristic of biomedical titanium alloy based on picosecond laser micromachining[J].Materials Science Forum，2018，939：104-109.[29] Schnell G，Staehlke S，Duenow U，et al. Femtosecond laser nano/micro textured Ti6Al4V surfaces—effect on wetting and MG-63 cell adhesion[J].Materials，2019，12（13）：2210.[30] Li C，Yang Y，Yang L J，et al. In vitro bioactivity and biocompatibility of bio-inspired Ti-6Al-4V alloy surfaces modified by combined laser micro/nano structuring[J].Molecules，2020，25（7）：1494.[31] Oberringer M，Akman E，Lee J，et al. Reduced myofibroblast differentiation on femtosecond laser treated 316LS stainless steel[J].Materials Science and Engineering：C，2013， 33（2）：901-908.[32] Shaikh S，Singh D，Subramanian M，et al. Femtosecond laser induced surface modification for prevention of bacterial adhesion on 45S5 bioactive glass[J].Journal of Non-Crystalline Solids，2018，482：63-72.[33] Yiannakou C，Simitzi C，Manousaki A，et al. Cell patterning via laser micro/nano structured silicon surfaces[J].Biofabrication，2017，9（2）：025024.[34] Michaljaničová I，Slepička P，Rimpelová S，et al. Regular pattern formation on surface of aromatic polymers and its cytocompatibility[J].Applied Surface Science，2016，370：131-141.[35] Wu H，Liu T，Xu Z Y，et al. Enhanced bacteriostatic activity，osteogenesis and osseointegration of silicon nitride/ polyetherketoneketone composites with femtosecond laser induced micro/nano structural surface[J].Applied Materials Today，2020，18：100523.（上接第84页）第30卷第4期2022年8月V ol.30 No.4Aug.2022安徽建筑大学学报Journal of Anhui Jianzhu UniversityDOI：10.11921/j.issn.2095-8382.20220417安徽宏村古村落非物质文化遗产数字化展示研究唐杰晓1，2，沈琛3（1.合肥师范学院艺术传媒学院，安徽合肥 230601；2.朝鲜大学设计大学院，韩国光州 61452；3.合肥师范学院计算机学院，安徽合肥 230601）摘‌要：安徽古村落非物质文化遗产是中华传统优秀文化的重要组成部分。