数学建模高斯扩散模型

§4-2高斯扩散模式ū —平均风速；Q—源强是指污染物排放速率;与空气中污染物质的浓度成正比,它是研究空气污染问题的基础数据;通常：ⅰ瞬时点源的源强以一次释放的总量表示；ⅱ连续点源以单位时间的释放量表示；ⅲ连续线源以单位时间单位长度的排放量表示；ⅳ连续面源以单位时间单位面积的排放量表示;δy—侧向扩散参数,污染物在y方向分布的标准偏差,是距离y的函数,m；δz—竖向扩散参数,污染物在z方向分布的标准偏差,是距离z的函数,m；未知量—浓度c、待定函数Ax、待定系数a、b；式①、②、③、④组成一方程组,四个方程式有四个未知数,故方程式可解;二、高斯扩散模式一连续点源的扩散连续点源一般指排放大量污染物的烟囱、放散管、通风口等;排放口安置在地面的称为地面点源,处于高空位置的称为高架点源;1. 大空间点源扩散高斯扩散公式的建立有如下假设：①风的平均流场稳定,风速均匀,风向平直；②污染物的浓度在y、z轴方向符合正态分布；③污染物在输送扩散中质量守恒；④污染源的源强均匀、连续;图5－9所示为点源的高斯扩散模式示意图;有效源位于坐标原点o处,平均风向与x 轴平行,并与x轴正向同向;假设点源在没有任何障碍物的自由空间扩散,不考虑下垫面的存在;大气中的扩散是具有y与z两个坐标方向的二维正态分布,当两坐标方向的随机变量独立时,分布密度为每个坐标方向的一维正态分布密度函数的乘积;由正态分布的假设条件②,参照正态分布函数的基本形式式5－15,取μ＝0,则在点源下风向任一点的浓度分布函数为：5－16式中 C—空间点x,y,z的污染物的浓度,mg/m3；Ax—待定函数；σy、σz—分别为水平、垂直方向的标准差,即y、x方向的扩散参数,m;由守恒和连续假设条件③和④,在任一垂直于x轴的烟流截面上有：5－17式中 q—源强,即单位时间内排放的污染物,μg/s；u—平均风速,m/s;将式5－16代入式5－17, 由风速稳定假设条件①,A与y、z无关,考虑到③和④,积分可得待定函数Ax：5－18将式5－18代入式5－16,得大空间连续点源的高斯扩散模式5－19式中,扩散系数σy、σz与大气稳定度和水平距离x有关,并随x的增大而增加;当y＝0,z ＝0时,Ax＝Cx,0,0,即Ax为x轴上的浓度,也是垂直于x轴截面上污染物的最大浓度点C max;当x→∞,σy及σz→∞,则C→0,表明污染物以在大气中得以完全扩散;2．高架点源扩散在点源的实际扩散中,污染物可能受到地面障碍物的阻挡,因此应当考虑地面对扩散的影响;处理的方法是,或者假定污染物在扩散过程中的质量不变,到达地面时不发生沉降或化学反应而全部反射；或者污染物在没有反射而被全部吸收,实际情况应在这两者之间;1高架点源扩散模式;点源在地面上的投影点o作为坐标原点,有效源位于z轴上某点, z＝H;高架有效源的高度由两部分组成,即H＝h＋Δh,其中h为排放口的有效高度,Δh是热烟流的浮升力和烟气以一定速度竖直离开排放口的冲力使烟流抬升的一个附加高度,如图5－10所示;当污染物到达地面后被全部反射时,可以按照全反射原理,用“像源法”来求解空间某点k的浓度;图5－10中k点的浓度显然比大空间点源扩散公式5－19计算值大,它是位于0,0,H的实源在k点扩散的浓度和反射回来的浓度的叠加;反射浓度可视为由一与实源对称的位于0,0,－H的像源假想源扩散到k点的浓度;由图可见,k点在以实源为原点的坐标系中的垂直坐标为z-H,则实源在k点扩散的浓度为式5－19的坐标沿z轴向下平移距离H：5－20k点在以像源为原点的坐标系中的垂直坐标为z＋H,则像源在k点扩散的浓度为式5－19的坐标沿z轴向上平移距离H：5－21由此,实源C s与像源C x之和即为k点的实际污染物浓度：5－22 若污染物到达地面后被完全吸收,则C x＝0,污染物浓度Cx,y,z,H＝C s,即式5－20; 2地面全部反射时的地面浓度;实际中,高架点源扩散问题中最关心的是地面浓度的分布状况,尤其是地面最大浓度值和它离源头的距离;在式5－22中,令z＝0,可得高架点源的地面浓度公式：5－23上式中进一步令y＝0则可得到沿x轴线上的浓度分布：5－24地面浓度分布如图图5－11所示;y方向的浓度以x轴为对称轴按正态分布；沿x轴线上,在污染物排放源附近地面浓度接近于零,然后顺风向不断增大,在离源一定距离时的某处,地面轴线上的浓度达到最大值,以后又逐渐减小;地面最大浓度值C max及其离源的距离x max可以由式5－24求导并取极值得到;令,由于σy、σz均为x的未知函数,最简单的情况可假定σy/σz＝常数,则当5－25 时,得地面浓度最大值5－26由式5－25可以看出,有效源H越高, x max处的σz值越大,而σz∝x max,则C max出现的位置离污染源的距离越远;式5－26表明,地面上最大浓度C max与有效源高度的平方及平均风速成反比,增加H可以有效地防止污染物在地面某一局部区域的聚积;式5－25和式5－26是在估算大气污染时经常选用的计算公式;由于它们是在σy/σz＝常数的假定下得到的,应用于小尺度湍流扩散更合适;除了极稳定或极不稳定的大气条件,通常可设σy/σz＝2估算最大地面浓度,其估算值与孤立高架点源如电厂烟囱附近的环境监测数据比较一致;通过理论或经验的方法可得σz＝fx的具体表达式,代入5－25可求出最大浓度点离源的距离x max,具体可查阅我国GB3840—91制定地方大气污染物排放标准的技术方法;3. 地面点源扩散对于地面点源,则有效源高度H＝0;当污染物到达地面后被全部反射时,可令式5－22中H＝0,即得出地面连续点源的高斯扩散公式:5－27其浓度是大空间连续点源扩散式5－19或地面无反射高架点源扩散式5－20在H＝0时的两倍,说明烟流的下半部分完全对称反射到上部分,使得浓度加倍;若取y与z等于零,则可得到沿x轴线上的浓度分布：5－28如果污染物到达地面后被完全吸收,其浓度即为地面无反射高架点源扩散式5－20在H＝0时的浓度,也即大空间连续点源扩散式5－19;高斯扩散模式的一般适用条件是：①地面开阔平坦,性质均匀,下垫面以上大气湍流稳定；②扩散处于同一大气温度层结中,扩散范围小于10km；③扩散物质随空气一起运动,在扩散输送过程中不产生化学反应,地面也不吸收污染物而全反射；④平均风向和风速平直稳定,且u＞1～2m/s;高斯扩散模式适应大气湍流的性质,物理概念明确,估算污染浓度的结果基本上能与实验资料相吻合,且只需利用常规气象资料即可进行简单的数学运算,因此使用最为普遍;二连续线源的扩散当污染物沿一水平方向连续排放时,可将其视为一线源,如汽车行驶在平坦开阔的公路上;线源在横风向排放的污染物浓度相等,这样,可将点源扩散的高斯模式对变量y积分,即可获得线源的高斯扩散模式;但由于线源排放路径相对固定,具有方向性,若取平均风向为x轴,则线源与平均风向未必同向;所以线源的情况较复杂,应当考虑线源与风向夹角以及线源的长度等问题;如果风向和线源的夹角β＞45,无限长连续线源下风向地面浓度分布为：5－29当β＜45时,以上模式不能应用;如果风向和线源的夹角垂直,即β＝90,可得：5－30对于有限长的线源,线源末端引起的“边缘效应”将对污染物的浓度分布有很大影响;随着污染物接受点距线源的距离增加,“边源效应”将在横风向距离的更远处起作用;因此在估算有限长污染源形成的浓度分布时,“边源效应”不能忽视;对于横风向的有限长线源,应以污染物接受点的平均风向为x轴;若线源的范围是从y1到y2,且y1＜y2,则有限长线源地面浓度分布为：5－31式中,s1＝y1/σy,s2＝y2/σy,积分值可从正态概率表中查出;三连续面源的扩散当众多的污染源在一地区内排放时,如城市中家庭炉灶的排放,可将它们作为面源来处理;因为这些污染源排放量很小但数量很大,若依点源来处理,将是非常繁杂的计算工作;常用的面源扩散模式为虚拟点源法,即将城市按污染源的分布和高低不同划分为若干个正方形,每一正方形视为一个面源单元,边长一般在～10km之间选取;这种方法假设：①有一距离为x0的虚拟点源位于面源单元形心的上风处,如图5-12所示,它在面源单元中心线处产生的烟流宽度为2y0＝σy0,等于面源单元宽度B；②面源单元向下风向扩散的浓度可用虚拟点源在下风向造成的同样的浓度所代替;根据污染物在面源范围内的分布状况,可分为以下两种虚拟点源扩散模式：第一种扩散模式假定污染物排放量集中在各面源单元的形心上;由假设①可得：5－32由确定的大气稳定度级别和上式求出的,应用P－G曲线图见下节可查取x o;再由x0＋x分布查出σy和σz,则面源下风向任一处的地面浓度由下式确定：5－33上式即为点源扩散的高斯模式5－24,式中H取面源的平均高度,m;如果排放源相对较高,而且高度相差较大,也可假定z方向上有一虚拟点源,由源的最初垂直分布的标准差确定,再由求出,由求出σz,由x0＋x 求出σy,最后代入式5－33求出地面浓度;第二种扩散模式假定污染物浓度均匀分布在面源的y方向,且扩散后的污染物全都均匀分布在长为πx0＋x／8的弧上,如图5-12所示;因此,利用式5－32求σy后,由稳定度级别应用P－G曲线图查出x0,再由x0＋x 查出σz,则面源下风向任一点的地面浓度由下式确定：5－34三、扩散参数及烟流抬升高度的确定高斯扩散公式的应用效果依赖于公式中的各个参数的准确程度,尤其是扩散参数σσz及烟流抬升高度Δh的估算;其中,平均风速u取多年观测的常规气象数据；源强q y、可以计算或测定,而σy、σz及Δh与气象条件和地面状况密切相关;1. 扩散参数σy、σz的估算扩散参数σy、σz是表示扩散范围及速率大小的特征量,也即正态分布函数的标准差;为了能较符合实际地确定这些扩散参数,许多研究工作致力于把浓度场和气象条件结合起来,提出了各种符合实验条件的扩散参数估计方法;其中应用较多的由是帕斯奎尔Pasquill 和吉福特Gifford提出的扩散参数估算方法,也称为P－G扩散曲线,如图5－13和图5－14所示;由图可见,只要利用当地常规气象观测资料,由表5－1查取帕斯奎尔大气稳定度等级,即可确定扩散参数;扩散参数σ具有如下规律：①σ随着离源距离增加而增大；②不稳定大气状态时的σ值大于稳定大气状态,因此大气湍流运动愈强,σ值愈大；③以上两种条件相同时,粗糙地面上的σ值大于平坦地面;由于利用常规气象资料便能确定帕斯奎尔大气稳定度,因此P－G扩散曲线简便实用;但是,P－G扩散曲线是利用观测资料统计结合理论分析得到的,其应用具有一定的经验性和局限性;σy是利用风向脉动资料和有限的扩散观测资料作出的推测估计,σz是在近距离应用了地面源在中性层结时的竖直扩散理论结果,也参照一些扩散试验资料后的推算,而稳定和强不稳定两种情况的数据纯系推测结果;一般,P－G扩散曲线较适用于近地源的小尺度扩散和开阔平坦的地形;实践表明,σy的近似估计与实际状况比较符合,但要对地面粗糙度和取样时间进行修正；σz的估计值与温度层结的关系很大,适用于近地源的lkm以内的扩散;因此,大气扩散参数的准确定量描述仍是深入研究的课题;估算地面最大浓度值C max及其离源的距离x max时,可先按式5－25计算出σz,并图5-14查取对应的x值,此值即为当时大气稳定度下的x max;然后从图5-13查取与x max对应的σy值,代如式5－26即可求出C max值;用该方法计算,在E、F级稳定度下误差较大,在D、C级时误差较小;H越高,误差越小;我国GB3840－91制定地方大气污染物排放标准的技术方法采用如下经验公式确定扩散参数σy、σz：5－35式中,γ1、α1、γ2及α2称为扩散系数;这些系数由实验确定,在一个相当长的x距离内为常数,可从GB3840－91的表中查取;。

大气污染扩散第一节大气结构与气象有效地防止大气污染的途径，除了采用除尘及废气净化装置等各种工程技术手段外，还需充分利用大气的湍流混合作用对污染物的扩散稀释能力，即大气的自净能力。

污染物从污染源排放到大气中的扩散过程及其危害程度，主要决定于气象因素，此外还与污染物的特征和排放特性，以及排放区的地形地貌状况有关。

下面简要介绍大气结构以及气象条件的一些基本概念。

一、大气的结构气象学中的大气是指地球引力作用下包围地球的空气层，其最外层的界限难以确定。

通常把自地面至1200 km左右范围内的空气层称做大气圈或大气层，而空气总质量的98．2％集中在距离地球表面30 km以下。

超过1200 km的范围，由于空气极其稀薄，一般视为宇宙空间。

自然状态的大气由多种气体的混合物、水蒸气和悬浮微粒组成。

其中，纯净干空气中的氧气、氮气和氩气三种主要成分的总和占空气体积的99.97％，它们之间的比例从地面直到90km高空基本不变，为大气的恒定的组分；二氧化碳由于燃料燃烧和动物的呼吸，陆地的含量比海上多，臭氧主要集中在55～60km高空，水蒸气含量在4％以下，在极地或沙漠区的体积分数接近于零，这些为大气的可变的组分；而来源于人类社会生产和火山爆发、森林火灾、海啸、地震等暂时性的灾害排放的煤烟、粉尘、氯化氢、硫化氢、硫氧化物、氮氧化物、碳氧化物为大气的不定的组分。

大气的结构是指垂直（即竖直）方向上大气的密度、温度及其组成的分布状况。

根据大气温度在垂直方向上的分布规律，可将大气划分为四层：对流层、平流层、中间层和暖层，如图5－1所示。

1. 对流层对流层是大气圈最靠近地面的一层，集中了大气质量的75％和几乎全部的水蒸气、微尘杂质。

受太阳辐射与大气环流的影响，对流层中空气的湍流运动和垂直方向混合比较强烈，主要的天气现象云雨风雪等都发生在这一层，有可能形成污染物易于扩散的气象条件，也可能生成对环境产生有危害的逆温气象条件。

因此，该层对大气污染物的扩散、输送和转化影响最大。

高斯扩散模型的适用条件1. 高斯扩散模型适用的条件之一就是要有相对稳定的环境呀！就好比在一个平静的湖泊里，水的流动很平稳，这时候高斯扩散模型就能很好地发挥作用啦！比如研究污染物在这样的环境中是怎么扩散的。

2. 它还适用于扩散源比较集中的情况呢！就像一个发光的灯泡，光线从那里散发出来，用高斯扩散模型来分析这种扩散是不是很合适呢？比如火灾中烟雾的扩散。

3. 扩散的物质不能有太奇怪的性质哦！可不是什么都能用高斯扩散模型的，这就像你不能用切菜的方法去绣花呀！比如一些特殊的化学物质可能就不太适用。

4. 要有足够的观测数据支持呀！没有数据就像巧妇难为无米之炊，怎么能让高斯扩散模型大展身手呢？比如对大气中颗粒物扩散的研究就得有大量数据。

5. 时间尺度也很重要呢！如果变化太快或太慢，高斯扩散模型可能就不太好使啦！好比一辆车开得太快或太慢，你都不好判断它的行驶轨迹，比如瞬间爆发的爆炸产生的扩散。

6. 空间范围也得合适呀！太大或太小的空间，高斯扩散模型也会有力不从心的时候呢！就像用小勺子舀大海的水，或者用大桶去装一滴水，比如研究小范围的气味扩散。

7. 系统不能太复杂啦！要是乱七八糟的因素太多，高斯扩散模型可就头疼咯！就像解一团乱麻，得先理清楚呀！比如生态系统中多种生物的相互作用下的物质扩散。

8. 扩散的速度得比较适中呀！太快或太慢，高斯扩散模型就不好把握啦！就像跑步，速度适中你才能更好地观察和分析，比如一些化学反应的扩散速度。

9. 环境不能总是变来变去的呀！一会儿这样一会儿那样，高斯扩散模型也会不知所措的！就像天气一会儿晴一会儿雨，怎么预测呀！比如海洋中水流和温度不断变化时的物质扩散。

10. 边界条件得明确呀！不然高斯扩散模型都不知道该从哪里开始从哪里结束呢！就像跑步没有起点和终点，怎么跑呀！比如研究一个房间内的气体扩散，房间的边界就得清楚。

我的观点结论就是：只有在这些条件满足的情况下，高斯扩散模型才能像一把锋利的宝剑，在研究扩散现象的战场上大显身手呀！。

高斯扩散模型 python高斯扩散模型是一种在空间或时间上描述随机现象扩散的模型。

该模型主要考虑了物质分子的运动规律，同时也考虑了物质分子之间的碰撞和相互作用。

在Python中，我们可以利用相关的库和函数来实现高斯扩散模型，下面将从以下几个方面进行介绍：1. 理解高斯分布首先，我们需要了解高斯分布，也称正态分布，表示连续变量的分布情况，常用于对连续变量进行建模和预测。

在Python中，可以使用SciPy库中的stats模块来计算高斯分布。

import numpy as npimport matplotlib.pyplot as pltfrom scipy.stats import norm# 定义均值(mean)和标准差(standard deviation)mu, sigma = 0, 0.1# 构造一些数据s = np.random.normal(mu, sigma, 1000)# 绘制直方图count, bins, ignored = plt.hist(s, 30, density=True)plt.plot(bins, 1/(sigma * np.sqrt(2 * np.pi)) *np.exp( - (bins - mu)**2 / (2 * sigma**2) ),linewidth=2, color='r')plt.show()2. 构建高斯扩散模型接下来，我们可以利用高斯分布的概念来构建高斯扩散模型。

在Python中，可以使用NumPy和SciPy库中的函数来进行计算。

import numpy as npimport matplotlib.pyplot as pltfrom scipy.special import erfcdef diffusion(x, t, D):return np.exp(-(x**2)/(4*D*t))/(2*np.sqrt(np.pi*D*t))# 初始化一些变量N = 100000 # 粒子数D = 1.0 # 扩散系数t = 100.0 # 时间bins = 200 # 直方图中的箱子数# 生成位移数据x = np.sqrt(2*D*t)*np.random.randn(N)# 利用扩散函数计算概率密度hist, bins = np.histogram(x, bins=bins, density=True) binCenters = 0.5*(bins[1:] + bins[:-1])theory = diffusion(binCenters, t, D)# 绘制直方图plt.plot(binCenters, hist, label='simulation')plt.plot(binCenters, theory, label='theory')plt.xlabel('x')plt.ylabel('P(x,t)')plt.legend()plt.show()3. 可视化高斯扩散模型的结果为了更直观地了解高斯扩散模型的结果，我们可以使用Matplotlib库中的plot函数来绘制直方图和散点图。

9.2.2大气污染物扩散的高斯模型模拟：可视化模拟点源大气污染的扩散9.2.2 Gaussian Atmospheric Dispersion Model突发性大气污染事故时有发生，对大气污染扩散进行模拟和分析，有利于减小事故的危害，减轻人员伤亡和财产损失。

高斯扩散模型是国际原子能机构(IAEA)推荐使用于重气云扩散模拟的数学模型,该模型在非重气云扩散的应用日益广泛。

高斯扩散模型是描述大气对有害气体的输移、扩散和稀释作用的物理或数学模型,是进行灾害预测和救援指挥的有力手段之一。

9.2.2.1高斯扩散模型高斯模型又分为高斯烟团模型和高斯烟羽模型。

大气污染物泄漏分为瞬时泄漏和连续泄漏，瞬时泄漏是指污染物泄放的时间相对于污染物扩散的时间较短如突发泄漏等的情形，连续泄漏则是指污染物泄放的时间较长的情形。

瞬时泄漏采用高斯烟团模型模拟，而连续泄漏采用高斯模型烟羽模型模拟。

高斯模型适用于非重气云气体，包括轻气云和中性气云气体。

要求气体在扩散过程中，风速均匀稳定。

在高斯烟团模型中，选择风向建立坐标系统，即取泄漏源为坐标原点，x 轴指向风向，y 轴表示在水平面内与风向垂直的方向，z 轴则指向与水平面垂直的方向，具体公式见式(9.1):22222222()()()22223/2(,,,)()(2)y x z z y x ut z H z H x y z Q C x y z t e e e e σσσσπσσσ--+----=⋅⋅⋅+⋅…………(9.1)其中：(,,,)C x y z t 为泄漏介质在某位置某时刻的浓度值；Q 为污染物单位时间排放量(mg/s); x σ、y σ、z σ分别x 、y 、z 轴上的扩散系数，需根据大气稳定度选择参数计算得到(m)；x 、y 、z 表示x 、y 、z 上的坐标值(m)；u 表示平均风速(m/s)；t 表示扩散时间(s)；H 表示泄漏源的高度(m)。

同理，高斯烟羽模型的表达式如：222222()()222(,,,)()2y z z y z H z H y z Q C x y z t e e e u σσσπσσ-+---=⋅⋅+………………………(9.2)9.2.2.2 技术方法若用高斯模型算出空间每一个点在一个时刻的污染浓度，这个计算量是很大的。

高斯扩散模型假设名词解释
高斯扩散模型是一种用来描述空气污染物在大气中传播和扩散
的数学模型。

它是基于高斯分布的假设，即空气污染物在水平方向上的传播服从正态分布。

在高斯扩散模型中，假设空气污染物在垂直方向上的传播是均匀的，即空气污染物在垂直方向上的浓度是恒定的。

这是基于大气中存在的湍流现象，使得空气混合均匀，污染物被均匀分散在大气中。

另外，高斯扩散模型还假设空气污染物在水平方向上的传播是径向对称的，即从污染源点开始，污染物浓度随着距离的增加呈现出高斯分布的特征。

这是因为在大气中存在着各种影响空气传播的因素，如风速、大气稳定度等，这些因素使得空气污染物向各个方向扩散。

高斯扩散模型可以通过一系列的数学公式来计算空气污染物在不同
位置的浓度分布。

这些公式考虑了污染源的排放强度、环境因素（如风速、大气稳定度等）以及地形特征等因素的影响。

通过模拟和计算，可以预测不同条件下空气污染物的传播范围和浓度分布，从而为环境管理和污染控制提供科学依据。

除了以上提到的假设，高斯扩散模型还可以考虑其他因素的影响，如地形地貌、建筑物的阻挡效应等，以更加准确地描述污染物在大气中
的传播过程。

它是环境科学领域中常用的一种模型，能够帮助我们更好地理解和管理空气污染问题。

第七节 扩散问题的偏微分方程模型物质的扩散问题，在石油开采、环境污染、疾病流行、化学反应、新闻传播、煤矿瓦斯爆炸、农田墒情、水利工程、生态问题、房屋基建、神经传导、药物在人体内分布以及超导、液晶、燃烧等诸多自然科学与工程技术领域，十分普遍地存在着. 显然，对这些问题的研究是十分必要的，其中的数学含量极大. 事实上，凡与反应扩散有关的现象，大都能由线性或非线性抛物型偏微分方程作为数学模型来定量或定性地加以解决.ＭＣＭ的试题来自实际，是“真问题⊕数学建模⊕计算机处理”的“三合一”准科研性质的一种竞赛，对上述这种有普遍意义和数学含量高，必须用计算机处理才能得到数值解的扩散问题，当然成为试题的重要来源，例如，ＡＭＣＭ-90Ａ，就是这类试题；ＡＭＣＭ-90Ａ要研究治疗帕金森症的多巴胺（dopamine ）在人脑中的分布，此药液注射后在脑子里经历的是扩散衰减过程，可以由线性抛物型方程这一数学模型来刻划. ＡＭＣＭ-90Ａ要研究单层住宅混凝土地板中的温度变化，也属扩散（热传导）问题，其数学模型与ＡＭＣＭ-90Ａ一样，也是线性抛物型方程.本文交代扩散问题建模的思路以及如何推导出相应的抛物型方程，如何利用积分变换求解、如何确定方程与解的表达式中的参数等关键数学过程，且以ＡＭＣＭ-90Ａ题为例，显示一个较细致的分析、建模、求解过程.§1 抛物型方程的导出设(,,,)u x y z t 是t 时刻点(,,)x y z 处一种物质的浓度. 任取一个闭曲面S ，它所围的区域是Ω，由于扩散，从t 到t t +∆时刻这段时间内，通过S 流入Ω的质量为2221(cos cos cos )dSd t ttSu u u M a b c t x y zαβγ+∆∂∂∂=++∂∂∂⎰⎰⎰. 由高斯公式得2222221222()d d d d t ttu u u M a b c x y z t x y z +∆Ω∂∂∂=++∂∂∂⎰⎰⎰⎰. （1） 其中，222,,a b c 分别是沿,,x y z 方向的扩散系数. 由于衰减（例如吸收、代谢等），Ω内的质量减少为22d d d d t ttM k u x y z t +∆Ω=⎰⎰⎰⎰，（2） 其中2k 是衰减系数.由物质不灭定律，在Ω内由于扩散与衰减的合作用，积存于Ω内的质量为12M M -.换一种角度看，Ω内由于深度之变化引起的质量增加为3[(,,,)(,,,)]d d d d d d d . (3)t ttM u x y z t t u x y z t x y zux y z t t Ω+∆Ω=+∆-∂=∂⎰⎰⎰⎰⎰⎰⎰显然312M M M =-，即2222222222d d d d ()d d d d .t ttt ttux y z t t u u u a b c k u x y z t x y z+∆Ω+∆Ω∂∂∂∂∂=++-∂∂∂⎰⎰⎰⎰⎰⎰⎰⎰由,,t t ∆Ω之任意性得2222222222u u u u a b c k u t x y z∂∂∂∂=++-∂∂∂∂ （4） 方程（４）是常系数线性抛物型方程，它就是有衰减的扩散过程的数学模型，对于具体问题，尚需与相应的定解条件（初始条件与边界条件等）匹配才能求得确定情况下的解.§2 Ｄirac 函数物理学家Dirac 为了物理模型之需要，硬是引入了一个当时颇遭微词的，使得数学与物理学传统密切关系出现裂痕的“怪”函数：0,0,() ()1.,0,x x x dx x δδ+∞-∞≠⎧==⎨∞=⎩⎰ （5）它的背景是清晰的，以一条无穷长的杆子为例，沿杆建立了一维坐标系，点的坐标为x ，杆的线密度是()x ρ，在(,]x -∞段，杆子质量为()m x ，则有d ()(), ()d ().d x m x x x x m x xρρ-∞==⎰. （6）设此无穷长的杆子总质量为１，质量集中在0x x =点，则应有001,,()0,,x x m x x x >⎧=⎨<⎩ 或写成 0()()m x H x x =-，其中()H x 为1,0,()0,0,x H x x >⎧=⎨<⎩ 如果沿用（6）中的算法，则在质量集中分布的这种情形有00,,(),0.x x x x ρ≠⎧=⎨∞=⎩且0()d ()xx x H x x ρ-∞=-⎰，于是得()d 1.x x ρ+∞-∞=⎰. （7）但是，从传统数学观点看，若一个函数除某点处处为零，则不论哪种意义下的积分，都必定为零，（7）式岂能成立！但是，δ函数对于物理学而言是如此之有用，以致物理学家正当地拒绝放弃它. 尽管当时数学家们大都嘲笑这种函数，但Ｐ.Ａ.Ｍ.Ｄirac 及其追随者们在物理领域却收获颇丰，Dirac 于1933年获诺贝尔物理奖. 当然Dirac 也意识到()x δ不是一个通常的函数，至于找一种什么办法来阐明()x δ这一符号的合法性，那就是数学家的任务了. 1940年，法国数学家许瓦兹（L.Schwartz ）严格证明了应用()x δ的正确性，把δ函数置于坚实的数学基础上；1950年，L. Schwartz 获数学界最高奖Fields 奖.δ函数的重要性质有：1）0()d 1x x x δ+∞-∞-=⎰. （8）2）00()()d ()x x f x x f x δ+∞-∞-=⎰. （9）其中()(,)f x C ∈-∞+∞，即0()x x δ-摘出了()f x 在0x x =的值.3）00()()dH x x x x dxδ-=-.（10） 4）()x δ的导数是存在的，不过要到积分号下去理解：00()()(),x x f x dx f x δ+∞-∞''-=-⎰ （11）()()00()()(1)().n n n x x f x dx f x δ+∞-∞-=-⎰（12）事实上，由于0()x x δ-在,+∞-∞处为零，则形式地用分部积分公式000()()()()d ()()d ,x x f x x x f x xx x f x x δδδ+∞+∞-∞-∞+∞-∞'---'=-⎰⎰其中，()(,)nf x C ∈-∞+∞，于是有（11）与（12）公式.5）对于()(,)x C ϕ∈-∞+∞，有000()()()()x x x x x x ϕδϕδ-=-.（13）6）1()() (0)||bx x b b δδ=≠.（14）7）000000(,,)()()()x x y y z z x x y y z z δδδδ---=---. （15）8）付立叶变换00[()].i x y y e λδ--= （16）[()] 1.x δ= （17）11221122[()()][()][()].C x x C x x C x x C x x δδδδ-+-=-+- （18）9）拉普拉斯变换00[(),[() 1.x x x e x δδδ--== （19）11221122[()]()][()[()].C x x C x x C x x C x x δδδδ-+-=-+- （20） 从上面的定义与性质看出，Delta 函数()x δ与一般可微函数还是有重大区别的，我们说它是“广义函数. ”§3 Cauchy 问题的解设扩散源在点000(,,)x y z 处，则此扩散问题满足Cauchy 问题2222222222000, (21)(,,,0)()()(). (22)u u u u a b c k u tx y z u x y z M x x y y z z δδδ⎧∂∂∂∂=++-⎪∂∂∂∂⎨⎪=---⎩对（21）（22）进行付立叶变换，且令123ˆ(,,), (,)[(,,,)]ut u x y z t λλλλλ==， 由于222222123222ˆˆˆ[], [], [],u u u uu u x y zλλλ∂∂∂=-=-=-∂∂∂ 102030000()[(,,,0)][()][()][()] ,i x y z u x y z M x x y y z z Me λλλδδδ-++=---= 故得常微分方程Cauchy 问题1020302222222123()ˆ()0,ˆ(0,).i x y z du a b c k udt u Meλλλλλλλ-++⎧++++=⎪⎨⎪=⎩ 得唯一解2222222123102030()()ˆ(,)a b c k t ix y z ut Me λλλλλλλ-+++-++=. （23）对（23）求逆变换1-，由于2122214[]a xa eλ---=， 2110221240[]()i x e aaex x λλ----=-， 故得12222000222ˆ(,,,)[]()()()exp 444u x y z t u x x y y z z k t a t b t c t -=⎧⎫---=----⎨⎬⎩⎭2222000222()()().444x x y y z z k t a t b t c t ⎧⎫---=----⎨⎬⎩⎭（24） 如果认为经过了相当长时间后，扩散已经终止，物质分布处于平衡状态，则方程（4）中的0ut∂=∂，于是有线性椭圆型方程的边值问题 22222222220, (,,)(,,)(,,).D u u u a b c k u x y z D xy z u x y z x y z ϕ∂⎧∂∂∂++-=∈⎪∂∂∂⎨⎪=⎩也可以用付立叶变换求解. 当然，根据实际情况，还可以考虑第二边条件(,,)Dux y z n ∂∂=ψ∂或第三边条件[](,,)D uu x y z nαβρ∂∂+=∂等，其中D ∂是区域D 的边界，n 是外法线方向，,αβ是实常数.§4 参数估计在Cauchy 问题（21）（22）的解（23）中，有四个未知的参数,,,a b c k ，它们分别是扩散与衰减过程中的扩散系数与衰减系数的算术平方根. 至于点源的质量与位置000,(,,)M x y z 是已知的.设观测取样为：11112222(,,,), (,,,),,(,,,),n n n n x y z m x y z m x y z m取样时刻为1t =（不然设00, t t t τ=是取样时间，则（21）变成2200t xx yy U t a U t b U =++2200zz t c U t k U -，对τ而言，取样时间为１,而方程形状与（21）一致），把在(,,)i i i x y z 点观测到的物质密度i m 与公式（24）都取对数，令1t =，则2222000222()()()ln (,,,1)ln []444x x y y z z u x y z abc k a b c ---=--+++. （25） 令222000222()()()111,,,,,,444x x y y z z X Y Z a b c αβγ---====-=-=-2ln ln abc k ε=--，则（25）写成 ln (,,,1)W u x y z X Y Z αβγε==+++，（26） 而我们已观测得(,,,)1,2,,i i i i X Y Z W i n =的数据，用三元回归分析方法求出,,,αβγε的估计值如下：ˆˆˆˆ()W X Y Z εαβγ=-++， （27）其中11111111, , , ,n n n nk i i i k k k k W W X X Y Y Z Z n n n n ========∑∑∑∑ˆˆˆ,,αβγ满足方程组 111213102122232031323330ˆˆˆ,,ˆˆˆ,,ˆˆˆ,.l l l l l l l l l l l l αβγαβγαβγ⎧++=⎪⎪++=⎨⎪++=⎪⎩其中10201130122211223311112131123211()(), ()(),()(),(), (), (),()(), ()(),()(), nnk k k k k k nk k k nn nk k k k k k nnk k k k k k nk k k l W X W W l Y Y W W l Z Z W W l X X l Y Y l Z Z l X X Y Y l X X Z Z l Y Y Z Z l ==========--=--=--=-=-=-=--=--=--∑∑∑∑∑∑∑∑∑1231133223, , .l l l l l ===由ˆˆˆ,,αβγ可求得222,,a b c 的估计值，即222111ˆˆˆ, , ˆa b cαβγ=-=-=-. 又由于2ln k abc ε=+- （28） 由（27）式可得ˆε，再把ˆˆˆ,,a b c 代入（28）得2ˆˆˆˆˆln kabc ε=+- （29）至此得到参数2222,,,a b c k 的估计值2222ˆˆˆˆ,,,a b c k ，把它们代入（24）分别替代2222,,,a b c k ，则得不含未知参数的解(,,,)u x y z t 的近似表达式.§5 竞赛试题分析AMCM-90A 不可用本文的思路与方法加以解决；该试题由东华盛顿大学数学系Ｙves Nievergelt 提供，要求研究药物在脑中的分布，题文称：“研究脑功能失调的人员欲测试新的药物的效果，例如治疗帕金森症往脑部注射多巴胺（Dopamine ）的效果，为了精确估计药物影响到的脑部区域，它们必须估计注射后药物在脑内空间分布区域的大小和形状.“研究数据包括50个圆柱体组织样本的每个样本药物含量的测定值（如图6-1），每个圆柱体长0.76mm ，直径0.66mm ，这些互相平行的圆柱体样本的中心位于网格距为1mm ×0.76×mm ×1mm 的格点上，所以圆柱体互相间在底面上接触，侧面互不接触. 注射是在最高计数的那个圆柱体的中心附近进行的. 自然在圆柱体之间以及由圆柱体样本的覆盖的区域外也有药物.“试估计受到药物影响的区域中药物的分布. ”“一个单位表示一个闪烁微粒的计数，或多巴胺的4.753×10-18克分子量，例如表6-1指出位于后排当中那个圆柱体的含药量是28353个单位. ”后方垂直截面164 442 1320 414 188 480 7022 14411 5158 352 2091 23027 28353 13138 681 7892126020921117317272131303 3765 1715 453前方垂直截面163 324 432 243 166 712 1055 6098 1048 232 2137 15531 19742 4785 330 444 11431 14960 3182 301 29420611036 258188图6-1数学模型只是实际问题的近似，要建立数学模型，一般首先要对所研究的实际问题进行必要和允许的简化与假设，而且，不同的简化与假设，又可能导致不同的数学模型，例如[2]是抛物型方程模型，而[3]则是椭圆方程模型.假设：（1）注射前大脑中的多巴胺含量可以忽略不计.（2）大脑中多巴胺注射液经历着扩散与衰减的过程，且沿,,x y z 三个方向的扩散系数分别是常数，衰减使质量之减少与深度成正比.（3）注射点在后排中央那个圆柱中心，即注射点的坐标000(,,)x y z 已知，注射量有医疗记录可查，是已知的.（4）注射瞬间完成，可视为点源delta 函数. （5）取样也是瞬间完成，取样时间已知为1t =.（6）样本区域与整个大脑相比可以忽略，样本组织远离脑之边界，不受大脑边界面的影响.在以上假设之下，显然可以用本文前面讲过的思路来建模，于是得AMCM-90A 的数学模型为Cauchy 问题（21）（22），解的表达式为（24），且用三元回归分析来估出参数,,,a b c k ，于是可以求得任意位置任意时刻药物的深度.如果所给数据认为是在平衡状态测得的，药物注射进脑后，从高深度处向低深度处扩散，与扩散同时，一部分药物进入脑细胞被吸收固定，扩散系数与吸收系数都是常数，但过一段时间，所有药物都被脑细胞所固定，达到了平衡态. 在这种假设下，[3]给出了下述的分析、建模、求解过程.设(,,,)v x y z t 是t 时刻在(,,)x y z 点处游离的药物浓度，(,,,)w x y z t 是t 时刻(,,)x y z 点处吸收固定的药物浓度，(,,)u x y z 是达到平衡态时(,,)x y z 点处吸收固定的药物浓度. 又设游离药物在各方向上有相同的扩散系数k ，吸收系数为h ，于是有vk v hv t∂=∆-∂. （30）又whv t∂=∂，即吸收速度与游离的浓度成正比，代入（30）得 ()v k w w t h t t∂∂∂=∆-∂∂∂. （31） 对（31）关于t 从０到+∞积分得t t t k vw wh+∞+∞+∞====∆-. （32）由于最后无游离药物，故(,,,)0v x y z +∞=，又开始时(0)t =无被吸收的药物，故(,,,0)0, (,,,0)0w x y z w x y z =∆=；平衡状态在t =+∞时达到，这时(,,)u x y z =(,,,)w x y z +∞，于是由（32）得(,,,0)ku u v x y z h-∆+=， （33）其中(,,,0)v x y z 是开始时的浓度分布，近似于注射点的点源脉冲函数. 把此注射点取为坐标原点(0,0,0)，则(,,,0)(,,),v x y z L x y z L δ=是注射量，于是2k h σ⎛⎫= ⎪⎝⎭记2(,,)u u L x y z σδ-∆+=， （34）作付立叶变换得22222222ˆˆ(),ˆ,1()s u u L Lus σξησξη+++==+++ 再作反变换得u σ-=-， （35）其中C 是可计算常数.如果考虑各向不同性，设,,x y z 方向上扩散系数分别为222,,a b c ，注射点在000(,,)x y z ，则222222000222()()()u u u a b c u L x x y y z z x y z δδδ⎛⎫∂∂∂-+++=--- ⎪∂∂∂⎝⎭， 于是解为(,,)u x y z =exp 1⎧⎪-⎨⎪⎩，（36） （36）中的D 可计算常数.用前面类似的方法可以进行参数估计.在建模过程中，点源函数的使用显然与实况有差别；尤其是认为扩散系数与吸收系数都是常数，对于人脑这种有复杂结构的区域，这种假设与实际不会完全符合；夜间与白天（睡与醒）对这些系数有无影响？脑中各点这些系数是否有变？除时间位置应考虑外，可能还与药液浓度有关. 如此看来，脑内药液分布的数学模型很可能不是常系数线性偏微分方程，而是函数系数的线性微分方程甚至是非线性偏微分方程. 这时，其解不再能用封闭公式来表达，求解过程会变得极为复杂，所以也可以考虑是否试用其他数学模型来解，例如在平衡态的假设下，用回归分析方法建立药液的模拟分布(,,)u f x y z =.对一个实际问题，其数学模型未必唯一，各模型间孰优孰劣，没有一般的判别法，须经实践来检验.参考文献[1]叶其孝，大学生数学建模竞赛辅导教材，湖南教育出版社，1993.[2]Christopher, R. Malone, Gian Pauletto, James, I. Zoellick, Distribution of Dopamine in the Brain, The Journal of Under graduate Mathematics, and its Applications, vol. 12(1991), Special Issue: The 1991 Mathematical Contest in Modeling, pp. 211-223.[3]孙晓东，荆秦，梁俊，脑中药物分布的数学模型，数学的实践与认识，1991年No. 4，63-69.[4]中国科学院数理统计组，常用数理统计方法，科学出版社，19784.。