精华内容
下载资源
问答
  • 用牛顿 同伦,参数微分法求非线性方程的根,其中里面有个子程序是牛顿法
  • 参数微分法中的欧拉法求非线性方程组的一组解
  • A.3 参数微分法

    2020-06-05 00:08:14
      对于系统x˙=f(t,x,p)(1)\dot{\mathbf{x}}=\mathbf{f}(t, \mathbf{x}, \mathbf{p})\tag{1}x˙=f(t,x,p)(1...是常数参数。   在许多情况下,式(1)的初值x(t0)x(t_0)x(t0​)以及p\mathbf{p}p都是不知道的,需要利用

      对于系统x˙=f(t,x,p)(1)\dot{\mathbf{x}}=\mathbf{f}(t, \mathbf{x}, \mathbf{p})\tag{1}

    式中p=[p1p2pq]T(2)\mathbf{p}=\left[p_{1} p_{2} \cdots p_{q}\right]^{T}\tag{2}

    是常数参数。
      在许多情况下,式(1)的初值x(t0)x(t_0)以及p\mathbf{p}都是不知道的,需要利用x(t)x(t)或其函数的测量值来进行估计。传统的估计方法(如最小二乘)需要求如下偏微分矩阵
    Φ(t,t0)=x(t)x(t0)(3)\Phi\left(t, t_{0}\right)=\frac{\partial \mathbf{x}(t)}{\partial \mathbf{x}\left(t_{0}\right)}\tag{3}

    Ψ(t,t0)=x(t)p\Psi\left(t, t_{0}\right)=\frac{\partial \mathbf{x}(t)}{\partial \mathbf{p}}
    因此,需要计算这两个微分矩阵。可获得这两个矩阵所应满足的微分方程
    Φ˙(t,t0)=F(t)Φ(t,t0),Φ(t0,t0)=IΨ˙(t,t0)=F(t)Ψ(t,t0)+f(t,x,p)p,Ψ(t0,t0)=0(4)\begin{array}{l} \qquad \dot{\Phi}\left(t, t_{0}\right)=F(t) \Phi\left(t, t_{0}\right), \quad \Phi\left(t_{0}, t_{0}\right)=I \\ \dot{\Psi}\left(t, t_{0}\right)=F(t) \Psi\left(t, t_{0}\right)+\frac{\partial \mathbf{f}(t, \mathbf{x}, \mathbf{p})}{\partial \mathbf{p}}, \quad \Psi\left(t_{0}, t_{0}\right)=0 \end{array}\tag{4}

    其中
    F(t)f(t,x,p)x(t)(5)F(t) \equiv \frac{\partial \mathbf{f}(t, \mathbf{x}, \mathbf{p})}{\partial \mathbf{x}(t)}\tag{5}

      微分矩阵存在如下关系
    δx(t)=Φ(t,t0)δx(t0)(6)\boldsymbol{\delta} \mathbf{x}(t)=\Phi\left(t, t_{0}\right) \boldsymbol{\delta} \mathbf{x}\left(t_{0}\right)\tag{6}

    其中δx\delta \mathbf{x}是偏离式(1)参考解xNx_N的一个小量。

    心得,评注

    1. 式(6)中的Φ(t,t0)\Phi\left(t, t_{0}\right)是对于参考解xNx_N轨迹而言的。当参考轨迹恒为0时,可得
      x(t)=Φ(t,t0)x(t0)(7) \mathbf{x}(t)=\Phi\left(t, t_{0}\right) \mathbf{x}\left(t_{0}\right)\tag{7}
    展开全文
  • 基于ENVI与ERDAS的Hyperion高光谱经验比值法、一阶微分法叶绿素及地表参数反演1 前期准备与本文理论部分1.1 几句闲谈1.2 背景知识1.2.1 Hyperion数据介绍1.2.2 遥感图像分类方法1.2.3 大气校正1.2.4 反演算法2 基于...

    1 前期准备与本文理论部分

    1.1 几句闲谈

      前面几篇博客介绍了基于Landsat这一多光谱遥感图像数据的多种地表温度(LST)反演方法,大家可以参考博客1博客2博客3;那么接下来,我们就将基于比多光谱数据可以说是更进一步的高光谱卫星数据——大名鼎鼎的Hyperion数据,进行多种其他地表参数的反演。其中,在此之前有必要先了解一下国内外主流的星载高光谱传感器及其平台的相关信息
      好啦,了解相关背景之后,话不多说,我们开始本次的内容。

      其中,本文所用遥感影像相关数据:
      点击链接:https://pan.baidu.com/s/16JUtCFhXlprRLjIhNBlikA
      提取码:csdn

    1.2 背景知识

    1.2.1 Hyperion数据介绍

      本文所使用的遥感影像数据不同于上述三篇博客所用的Landsat-7 ETM+影像数据,其来自地球观测卫星-1(Earth Observing-1,EO-1)。EO-1是美国国家航空航天局(National Aeronautics and Space Administration,NASA)新千年计划(New Millennium Program,NMP)地球探测部分中第一颗对地观测卫星,其目的即为在21世纪接替Landsat-7卫星,于2000年11月发射升空。除此之外,NMP目前还包括深空探测(Deep Space,DS)、空间技术(Space Technology,ST)两个部分。
      EO-1卫星轨道参数与Landsat-7较为近似,以期实现两颗卫星图像每天具有1至4景的重叠,从而进行二者的对比。此外,EO-1已于2017年2月停止服役。
      EO-1搭载了三种传感器,分别为高光谱成像光谱仪(Hyperion)、高级陆地成像仪(Advanced Land Imager,ALI)与线性标准成像光谱仪阵列大气校正器(the Linear Etalon Imaging Spectrometer Array Atmospheric Corrector,LAC)。一般地,传统的陆地资源卫星只能提供为数不多的多光谱波段,并不能很好满足日常实际研究、运用的需要;而借助具有242个波段、光谱范围为356至2578纳米的EO-1 Hyperion传感器,可获得更具价值的高光谱数据。
      EO-1遥感影像命名格式如下:
      EO1SPPPRRRYYYYDDDXXXML_GGG_VV
      其中,EO1为EO-1卫星代号,S为所用传感器代号(H为Hyperion传感器,A为ALI传感器),PPP为成像时目标所处全球参考系统(Worldwide Reference System,WRS)的轨道(Path),RRR为成像时目标所处WRS的行(Row),YYYY为成像年份,DDD为成像当日在该年份的天数,XXX分别为Hyperion、ALI与AC三种传感器的开关状态(1为开启,0为关闭),M为指向模式【(N为天底模式(Nadir),P为点在轨道模式(Pointed Within Path/Row),K为点不在轨道模式(Pointed Outside Path/Row)】,L为图像长度【F为全景(Full Scene),P为局部景(Partial Scene),Q为次级局部景(Second Partial Scene),S为样例(Swatch)】,GGG为影像地面接收站,VV为影像版本编号。
      一般地,遥感卫星传感器主要有两大类型:摆扫式(Whisk Broom Scanners)与推扫式(Push Broom Scanners);Hyperion属于后者。其242个波段分为可见光近红外波段(V-NIR)与短波红外波段(SWIR);其中,1至70波段属于V-NIR通道,71至242波段属于SWIR通道。两个波段之间具有20个波段的波长数值相互重叠,其分别用两套不同的敏感元件收集各自信号。
    Hyperion产品波段信息如下所示。
    在这里插入图片描述
      一般地,辐射校正包括辐射定标、大气校正和太阳及地形校正,用来消除辐射误差;而上述“辐射校正”包括正射校正,即使用地形数据的几何校正,不包括大气校正。
      Hyperion产品分为两级,Level 0与Level 1。前者为原始数据,其仅用来生产Level 1产品。Level 1产品则可以继续分为L1A、L1B、L1R、L1Gs与L1Gst等。其中,L1B产品与L1R产品分别由美国TRW与USGS处理生成。上述两种产品与L1A产品的最大不同在于,前二者纠正了V-NIR波段与SWIR波段的空间错位问题。
      而通过图像原始数据中“README.txt”文件可知,本文所使用的数据为L1Gst级别(大家自己操作时,也依据自己文件的实际情况判断即可);如以上内容所述,其已经过部分预处理过程,如正射校正、坐标转换等。
      Hyperion产品图像数据的空间分辨率为30米。
      如上述内容所提到的,Hyperion产品共有242个波段,波长覆盖范围为356至2578纳米。而在这其中,由于紫色光与短波红外两端有些通道的响应度较低,数据利用价值不大,因此系统并未对这些波段加以辐射校正——即对于正式分发的Hyperion产品Level 1数据,上述242个波段中,1至7波段、58至76波段、225至242波段亦称之坏波段(Bad Bands),其数值均为0值。除此之外,121至126(或127)波段、167至178(或180)波段、224(或222至224)波段由于受到水汽吸收影响较为严重,亦需要作为坏波段处理。最后,通过上述筛选后剩余的波段中,还包括波长重叠的56至57波段和77至78波段;考虑到77至78波段噪声相对前者较大,信噪比较低,因此将其同样剔除。
      上述波段筛选结果如下所示。
    在这里插入图片描述
      此外,除上述因未定标和水汽吸收影响被剔除的波段外,还包括坏线、条纹与“Smile”效应等非正常像元。其中,坏线是指一行或一列DN值为零或非常小的像元;条纹是指像元DN值不为零但较小,与周围有明显差异的带状现象;“Smile”效应是指由于前期光谱定标而产生的光谱差异。若对定量遥感结果精度要求较高,上述这些非正常像元或现象均需要在预处理阶段加以解决。
      目前,针对Hyperion产品图像的坏线处理主要通过均值邻域法,即用坏线周围的像元像素值求平均后替换坏线中的像素——这依赖于坏线像元与周围像元的高度相关关系;若当坏线出现在地形比较复杂的区域,采用这种方法则会导致处理结果精度较低。
      对此,有国内学者提出了一种基于地物类型的坏线修复方法。该方法假设各波段图像成像在一瞬间完成,或成像过程中天气未发生变化,则不同波段图像中同一类地物的DN值应当相同。当某一波段的图像A的地物M出现部分坏线时,搜索另一对应位置无坏线的波段对应的图像B,确定B中M位置对应的DN值α;随后继续在B中搜索,找到所有DN值为α的像元N。随后返回A,求取A中N像元对应位置的DN值,将其求平均后作为M位置处的DN值。这一坏线修复方法大大缩小了修复结果与真实值之间的均方根误差,效果较好。
      条纹尤多出现于SWIR波段,其严重影响图像质量;因此,条纹去除同样十分重要。针对这一问题,有学者提出了一种全局去除条纹的方法。其原理公式如下:
    在这里插入图片描述
      其中,〖x^’〗_ijk为第k波段对应图像第i行,第j列像素校正后的数值,x_ijk为其原值;s_ik为第k波段对应图像第i列标准差,m_ik为第k波段对应图像第i列平均值,(s_ik ) ̅与(m_ik ) ̅分别为第k波段对应图像第i列标准差与平均值的对应平均值。实际计算时,用整幅图像的平均值与标准差代替每一列的平均值与标准差。
      针对条纹去除结果的检验,可以利用最低噪声分离旋转(Minimum Noise Fraction Rotation,MNF)转换加以实现。
      “Smile”效应是指在垂直飞行方向上,像元的波长随着图像中心位置向两旁移动,产生的偏移。其中,V-NIR波段的“Smile”效应较为明显,约在2.6至3.6纳米范围内;SWIR波段的这一效应表现则不如前者明显,其约在0.40至0.97纳米范围内。“Smile”效应同样可以借助光谱信息的MNF转换加以检测。“Smile”效应往往会对植被遥感产生影响;可以通过大气纠正来降低“Smile”效应的影响。
      Hyperion产品Level 1数据以有符号的整型数据格式记录其信息,数值范围为-32768至32767;而实际地物辐射数值都比较小,如V-NIR波段最大辐射约为750 W/m2/sr/Lm,而SWIR波段最大辐射仅有350 W/m2/sr/Lm。因此,在产品生成时,系统采用了扩大因子,即分别将V-NIR与SWIR所对应波段的数值扩大为40和80倍。出于这个原因,在应用Hyperion产品Level 1数据时,需要将像元值转换为绝对辐射值。其中,两通道对应波段的计算公式有所不同:
    在这里插入图片描述
      其中,L_ℷVNIR为V-NIR波段定标后的绝对辐射值,L_ℷSWIR为SWIR波段定标后的绝对辐射值。
      此外,在计算时需要注意,由于经过波段筛选后的图像波段不再完全连续,会出现一些间段区域。因此需要首先将同一通道内波段编号相连的波段合并,再将同一通道内全部波段合成;第二次合成后,将两个通道对应的波段分别带入上述两个公式,计算结束后再将两个通道的176个波段最终合成为一幅图像。
      编辑头文件时,除需注意中心波长的修正外,还要注意半高宽(又称半峰全宽,Full Width at Half Maxima,FWHM)的输入。若不进行大气校正或进行大气校正但不采取FLAASH大气校正方法,则无需导入FWHM这一参数。
      需要注意的是,本文重点在于Hyperion的反演操作,因此上述Hyperion图像的坏线修复、条纹去除等处理暂未写入本文。具体操作我将在后期的博客中展现。

    1.2.2 遥感图像分类方法

      本文使用监督分类(又称训练分类,Supervised Classification)方法提取影像中各个不同地物部分,尤其是太湖湖面水体区域。而和前期博客一致,本文的监督分类部分依然借助ERDAS软件,选择采取平行六面体规则(Parallelepied)这一非参数规则(Non-parametric Rule)与最大似然规则(Maximum Likelihood)这一参数规则(Parametric Rule)加以实现;并且尝试了特征空间规则(Feature Space)的效果。
      监督分类是指在已掌握有足够先验知识(亦即训练场地)的情况下,根据已有训练场地提供的已知属性样本选择特征参数,并训练、建立得到对应判别函数;随后进而将图像未知类别部分的像素的值代入建立得到的判别函数,依据所选择的不同判别准则,对该样本所属的地物类别进行自动分类处理;即监督分类是一种利用已知地物属性、信息,进而对未知属性地物进行分类的方法。
      正如上述内容所提到的,常用的监督分类判别规则包括非参数规则与参数规则两个部分。其中,非参数规则包括特征空间规则与平行六面体规则;参数规则则包括最小距离规则(Minimum Distance)、马氏距离规则(Mahalanobis Distance)、最大似然规则以及波普角规则(Spectral Angle Mapper)等多种。
      ERDAS IMAGINE 2015软件中“Supervised Classification”模块可以选择的规则十分齐全。本次我们依然运用这一模块,并选择采取平行六面体规则与最大似然规则加以实现。平行六面体规则是指:根据训练样本的图像亮度值,形成一个N维的平行六面体数据空间;其余像元的光谱数值如果落在平行六面体任何一个训练样本所对应的区域,其就被划分至这一对应的类别中。其中,平行六面体的尺度是由标准差阈值所确定的,而该标准差阈值则是根据所选类的均值求出。最大似然规则是指:假设每一个波段的每一类统计都呈正态分布,并计算给定像元属于某一训练样本的似然度;随后,像元将最终被归并到似然度最大的一类样本当中。
      另一方面,遥感影像处理也可使用非监督分类(又称聚类分析或点群分类,Unsupervised Classification)加以完成。非监督分类方法是在多光谱图像中搜寻、定义其自然相似光谱集群的过程,其不必针对影像地物获取先验知识,仅依靠影像中不同类地物的光谱(或纹理)信息进行特征提取,再统计特征的差别,从而达到分类的目的;最后对已分出各个类别的实际属性进行确认。

    1.2.3 大气校正

      本次实验可以不进行大气校正;但不经过大气校正,由于未消除大气影响,会使得得到的结果会有一定误差。本实验中,我们首先不进行大气校正,对叶绿素加以反演;随后再先后尝试FLAASH大气校正与QUAC快速大气校正两种方法,对得到的高光谱遥感影像加以处理,从而对叶绿素反演精度加以提升。
      大气层是介于卫星传感器与地球表层之间的一层由多种气体及气溶胶组成的介质层,其由下至上可依次分为:对流层(0至8-18千米)、平流层(8-18至50-55千米)、中间层(50-55至80-85千米)、电离层(又称热层,80-85至800千米)、散逸层(800至3000千米)。在太阳辐射到达地表再到达卫星传感器这一过程中,辐射两次经过大气,故大气对太阳辐射的作用影响较大。
      大气校正的目的即消除大气和光照等因素对地物反射的影响,广义上讲是获得地物反射率、辐射率或地表亮度温度等真实物理模型参数,狭义上讲是获取地物真实反射率数据。大气校正可以用来消除大气中水蒸气、氧气、二氧化碳、甲烷和臭氧等物质对地物反射的影响,也可以消除大气分子和气溶胶散射的影响。在大多数情况下,大气校正也是反演地物真实反射率的过程。
      有时可完全忽略遥感数据的大气影响。对某些地物分类和变化检测过程(如用最大似然法对单时相遥感数据进行分类)而言,只要影像中用于分类的训练数据具有相对一致的尺度(即均经过或未经过大气校正),大气校正并不是必需的;而当训练数据需要进行时空拓展时,其影像分类和各种变化检测则需进行大气校正。
      而对于涉及定量遥感的应用方面,如前期博客所进行的地表真实温度反演运算,大气校正则多是必不可少的环节。有研究指出,是否进行大气校正过程,对植被稀少或植被被破坏地区的NDVI计算误差影响可达50%。
      大气校正方法整体上可以分为两类:一是绝对大气校正方法,即将遥感图像的DN值转换为地表反射率、地表辐射率、地表温度等数值的方法;二是相对大气校正方法,即将原始的DN值校正为不考虑地物实际反射率情况下,相同DN值代表相同地表反射率的格式的结果。在此其中,常用的绝对大气校正方法包括基于辐射传输模型方法(其中包括MOR/DTRAN模型、LOWTRAN模型、ATCOR模型、5S模型、6S模型等)、基于简化辐射传输模型的黑暗像元方法、基于统计学模型的反射率反演方法等;常用的相对大气校正方法包括基于统计的不变目标方法、直方图匹配方法等。针对上述繁多方法的选择,当进行精细的定量遥感研究时,需要采用基于辐射传输模型的大气校正方法;当进行动态监测研究时,可采用相对大气校正或简单的绝对大气校正方法;但所知参数较少时,则只能选择相对简单、对参数要求较少的方法。
      本文采用FLAASH大气校正与QUAC快速大气校正两种方法。
      FLAASH大气校正是一种绝对大气校正方法,其适用于多光谱数据与高光谱卫星影像数据,能够精确补偿大气对辐射的影响。这一方法采用 MODTRAN4+ 辐射传输模型,需要输入地表反射率,通过影像像素光谱中的特征估计大气属性,可以有效去除水蒸气、气溶胶散射等效应带来的干扰,精度较高。而另一方面,通过博客1可知,FLAASH大气校正方法同样对待处理数据有着一定严格的要求,例如需要数据的存储结构为“BIP”或“BIL”模式,像元值类型为经过定标后的辐射亮度(即辐射率)数据,数据单位为(μW)/(cm2nmsr)等。这需要我们针对原始数据加以一定处理。
      QUAC快速大气校正自动从图像中收集不同物质的波谱信息,获取经验值,从而完成高光谱和多光谱的快速大气校正。这一大气校正方法只可以对多光谱、高光谱图像进行处理。
      值得一提的是,FLAASH大气校正与QUAC快速大气校正两种方法得到的结果均为扩大了10000倍的反射率数据。

    1.2.4 反演算法

      获得太湖水体区域多个样点的高光谱数据曲线后,分别对其加以一阶微分处理与倒数对数处理。
      本文中,针对太湖水体叶绿素a含量的估计,首先采用常见文献中已给出的模型与参数:
    在这里插入图片描述
      其中,Chl.a为水体叶绿素a含量,单位为(mg/L);R_1与R_2为两个特征波段,R_1为近红外波段的反射峰,R_2为红光波段的吸收谷,二者需要通过光谱曲线加以求取;0.325与0.311为两个经验模型参数。一般地,特征光谱曲线其可以依据实测水体叶绿素a含量与光谱数值,通过计算二者相关系数随着对应波长变化而变化的关系,从而加以确定;而经验参数则是依据选择出的特征波段,通过回归分析等拟合手段求出。这一建立模型的正演过程本次我们不做讨论。
      另一方面,由于所获取的高光谱数据可能会受到背景噪声影像,且相似光谱之间因存在明显共线性而导致信息冗余,我们还可以对光谱数据加以一定预处理,从而消除上述噪声,突出光谱特征。往往使用光谱一阶微分、光谱倒数对数与光谱包络线去除等方式加以处理。其中,光谱包络线去除可以方便地在ENVI软件中实现。我们本次采取光谱一阶微分、光谱倒数对数方法。
      光谱一阶微分可以去除部分线性或接近线性的背景,以及噪声光谱对目标光谱的影响。光谱一阶微分计算公式如下:
    在这里插入图片描述
      其中,R^’ (λ_i )为波段λ_i对应的一阶微分数值,R(λ_i )为波段λ_i原本数值,λ_(i+1)为对应波段中心波长。
      在光谱定量分析中,光谱倒数对数可以十分有效地增强相似光谱之间的差别;计算如下:
    在这里插入图片描述
      其中,〖R(λ_i )〗_L为波段λ_i对应的倒数对数数值,R(λ_i )为波段λ_i原本数值,λ_i为对应波段中心波长。
      此外,包络线去除法亦是一种常用的光谱分析方法,其可以有效地突出光谱曲线的吸收和反射特征。但我们在此不做讨论,后期我会专门写一篇博客讨论这一内容。
      除水体叶绿素a含量反演外,利用高光谱数据反演地表沉积物颗粒度参数同样被学者们加以广泛研究。有学者借助Hyperion遥感影像,分别建立一个单波段因子与四个多波段组合因子;利用光谱包络线去除方法,挑选出与颗粒度参数相关性较高的几个特征波长,并选用稳定性更高、整体更简单的线性回归模型作为最终的参数反演模型。以平均粒径这一参数为例,具体反演计算模型如下:
    在这里插入图片描述
      其中,d_m为平均粒径,F为特征波段的运算函数,λ_1与λ_2分别为所选择的两个特征波段。其中,
    在这里插入图片描述
      同时,本次土壤有机质含量反演模型采用如下公式:
    在这里插入图片描述
      其中,OM为土壤有机质含量。

    2 基于经验比值法、一阶微分法的叶绿素a含量反演

      首先,在未进行大气校正的情况下,我们基于经验比值法、一阶微分法,对叶绿素a含量加以反演。

    2.1 数据导入与波段合成

    (1) 将数据中“EO1H1190382004095110KZ.tgz”图像压缩包文件解压缩;打开ENVI Classic 5.3(64-bit)软件,选择“File”→“Open Image File”,在弹出的文件选择窗口中分别选中压缩包中后缀名包括8至57的50个“.tif”格式文件;点击“打开”。

    在这里插入图片描述

    在这里插入图片描述
    (2) 选择→“Layer Stacking”,在弹出的文件选择窗口中选择打开的50个波段图像。

    在这里插入图片描述
    (3) 选择后,需要观察六幅图像在待合成文件列表中的排列顺序。由于后期我们需要对合成后的图像各个波段的信息进行折线图形式的分析操作,因此需要将待合成文件按照“中心波长值”由小至大的顺序排列。点击“Reorder Files”即可实现这一功能。

    在这里插入图片描述在这里插入图片描述

    (4) 结合本文开头提到的三篇博客可知,由于导入“Layer Stacking”模块的波段往往是倒序排列,这使得几十个波段的合成排序操作变得较为复杂——需要进行排序操作百余次。因此本次尝试再将文件导入ENVI软件时就采取倒序导入的方式。但随后发现即使这样操作,也会使得顺序列表中的波段倒序排列。
    在这里插入图片描述
    在这里插入图片描述

    (5) 顺序排列完毕后,检查投影信息等无误后,点击“OK”即可开始合成操作。
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

    (6) 多次重复这一过程,直至将同一通道内全部波长编号相连的图像各自拼接。
    在这里插入图片描述在这里插入图片描述
    在这里插入图片描述

    (7) 四个波合成后,将SWIR通道对应的三个一次合成后的图像再一次合成。
    在这里插入图片描述

    2.2 辐射定标与波段合成

    (1) 选择“Basic Tools”→“Spectral Math”,在弹出的公式创建窗口中输入本次实验的两个辐射定标公式;输入单个公式完成后,点击下方“Add to List”按钮,即可将键入的公式存入待选择区内。为了减少后期不必要的工作量,可以在编辑完成每一公式后点击“Save”按钮,并将最近一次保存的同一公式文件覆盖,以将待选择区内的两条公式统一保存。若直接选择“Band Math”(如下图)会导致无法对各波段实现简单的批量操作。

    在这里插入图片描述
    (2) 保存公式完成后,点击左下角“OK”按钮,即可开始公式的计算。在弹出的公式变量文件选择窗口中,将这一公式的变量“S1”选择为我们通过上述操作生成、添加的图像文件“layerstacking1_8_57”,并将第二个公式的变量“S1”选择为同样通过上述步骤获得、添加的图像文件“layerstacking5_79_223”。随后,配置对应两个定标输出文件地址等信息。
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述在这里插入图片描述

    (3) 配置完成后,点击左下角“OK”按钮,即可开始公式的运行。运行结束,将所得到的研究区两幅辐射定标结果图像导入ENVI软件中,显示如下。

    在这里插入图片描述
    (4) 由两幅图各自的光谱曲线可以看出,V-NIR波段数约为50个,SWIR波段数约为120个。这与我们通过上述波段合成得到的波段数分别为50个、126个的两幅图像相符合。
    (5) 再次执行上述波段合成操作,将V-NIR波段图像与SWIR波段图像合成。

    在这里插入图片描述
    (6) 通过本次波段合成,可以看到得到的结果图像包含了V-NIR波段数50个,SWIR波段数126个。

    2.3 编辑头文件

    (1) 由于在选择波段合成范围时,我们去除了一些“坏波段”,因此合成后的图像中心波长、FWHM等信息均丢失,需要重新将其导入。考虑到176个波段数量较大,我们可以采取Excel数据导出的方式将其填入ENVI影像中。
    (2) 首先在Microsoft Office Excel软件中打开原始遥感影像数据文件夹中的“Hyperion_Spectral_Coverage.xls”文件。这一表格文件为我们列出了Hyperion产品Level 1数据242个波段的相关信息,包括波段编号(Hyperion Band)、平均中心波长(Average Wavelength)、FWHM等属性数据。
    (3) 将我们最终定标、合成后的176个波段对应的上述三种属性数据导入到一张新的表格中,并导出为TXT文件。其中需要注意,将V-NIR波段与SWIR波段重叠部分的数据剔除。

    在这里插入图片描述
    (4) 其次在“Available Bands List”文件列表中右键选择波段合成后的结果文件,选择“Edit Header”功能。
    (5) 选择“Edit Attributes”→“Wavelengths”,将TXT文件导入,并注意将数据单位修改为“Micrometers”。
    在这里插入图片描述
    在这里插入图片描述

    (6) 此时再次打开上述操作编辑头文件完毕的176个波段的图像,可以看到光谱曲线中横坐标已经由原先的1至176转变为对应波段的中心波长数值范围。
    在这里插入图片描述

    2.4 图像格式转换

    经过ENVI处理后,得到的结果图像并不能直接被ERDAS软件所识别,我们需要将其转换为“.img”等格式文件。
    (1) 选择“File”→“Save File As”→“ERDAS IMAGINE”,在弹出的文件转换选择窗口中选择经过波段合成与头文件编辑后的结果图像,点击左下角“OK”。

    在这里插入图片描述
    (2) 在弹出的转换文件属性配置窗口中设置,配置好结果图像文件保存路径、保存文件名等。

    在这里插入图片描述
    (3) 文件格式转换完毕后,即可将结果文件在ERDAS软件中打开,并在其中进行后续操作。
    在这里插入图片描述

    2.5 EDRDAS文件导入与裁剪

    通过上述步骤得到的“.img”格式文件范围为太湖区域加之其西南部大面积陆地区域,水体占比并不多;而遥感图像整体区域面积较大,文件数据量多;另一方面,我们的主要反演目标区域为太湖水体。因此,我们需要借助自行划定的AOI(Area of Interest)文件,将原本的整体面积地区裁剪为太湖占比较多的地区。
    (1) 打开ERDAS IMAGINE 2015软件,在黑色图层窗口区域右键,并选择“Open Raster Layer”,在弹出的文件打开选择窗口中选择经过上述全部预处理后的“.img”结果图像,点击右上角“OK”。ERDAS软件选择文件时,可以借助“Recent”和“Goto”功能,选择最近操作的文件夹路径或文件。
    (2) 在图像中划定一个合适的太湖AOI区域,使得全部太湖水体都包含在内。
    (3) 选择“Raster”→“Subset & Chip”→“Create Subset Image”,在弹出的文件裁剪配置窗口中选择上述“.img”结果图像,配置输出文件路径和文件名,并选择输出图像文件格式为“Float Single”,在下方“AOI”选项中选择“Viewer”,点击左下角“OK”选项。这样即可以刚刚划定的AOI区域为裁剪范围。
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

    (4) 观察得到的结果图像,若对其直方图加以分析可发现,图像中存在大量数值为零的像元。其实这样并不是图像出现错误,而是在倾斜长方形图像周围区域实际上仍为图像所包括的区域,这些大量黑色的像素数值为0,使得直方图出现零值。
    在这里插入图片描述
    在这里插入图片描述

    2.6 监督分类

    如本文第一部分原理所述,本次实验我们依然利用监督分类方式区分太湖水体、其他水体与其它土地,从而实现更加精确的太湖区域提取。监督分类依然在ERDAS中实现。
    (1) 选择“Raster”→“Supervised”→“Supervised Editor”,在弹出的AOI区域显示表中可以看到,此时还没有添加进入任何AOI,表格中处于空白状态。
    (2) 分别依据太湖水体、其他水体与其它土地,在图像中划定区域。每一种地表类型需要划出尽可能多的AOI,以期提高最终监督分类的效果与精确度。
    在这里插入图片描述

    (3) 每划定一个AOI,便添加进入“Supervised Editor”中;同一地物类型的AOI添加完毕后,对这一类型的全部AOI加以合并,并赋以一个易于辨认的Value值。
    在这里插入图片描述
    在这里插入图片描述

    (4) 其中,可以借助卫星地图影像对具体地物类别加以具体区分。但应注意,我们所使用的遥感影像为2004年的数据,而卫星影像数据成像时间相对较为接近我们现在的时间。因此,利用卫星影像亦只可以作为一种参考。
    (5) 另一方面,亦发现无论是采用灰度图像还是选择三个波段分别作为R、G、B的值,在“Supervised Editor”中显示的颜色都是一致的,即原有的灰度图像。这一特点与前几次的多光谱博客有所差距,或许为高光谱影像的特征之一。

    在这里插入图片描述

    在这里插入图片描述
    (6) AOI划分完毕后,保存。选择“Raster”→“Supervised”→“Supervised Classification”,在弹出的监督分类配置窗口中配置输入文件、输出文件与用于指定区域地表类型的AOI文件。

    在这里插入图片描述
    (7) 得到的监督分类结果如下。其中,选择了采取平行六面体规则这一非参数规则。可以看到,由于三种地类分类颜色相近,具体的分类结果需要放大后才可以看清。
    在这里插入图片描述
    在这里插入图片描述

    (8) 下图中,左图则是采用特征空间规则这一非参数规则加以实现的监督分类结果。二者对比可以看出,右侧的平行六面体规则分类效果明显优于左侧。
    在这里插入图片描述

    2.7 水体光谱曲线提取

    (1) 利用“Inquire”模块,分别选取太湖水体中20处不同位置,并上述位置对应的各波段光谱数值导入Excel表格中,作为我们观察特征波段的基础。
    (2) 在进行这一步骤时,需要注意在ERDAS软件左下角有一个百分比显示区域。由于176个波段数量较多,数据在复制过程中会有些卡顿;借助百分比显示状态即可看出复制进度。
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述

    2.8 特征波段选取与计算

    (1) 将20个采样点的光谱曲线导入Excel软件后,利用本文第一部分原理中的公式计算对应光谱数值的一阶微分与倒数对数。并将得到的结果制图。以下三幅图分别为光谱原始数值、一阶微分与倒数对数对应图像。
    在这里插入图片描述
    在这里插入图片描述

    (2) 依据以上三幅光谱曲线,我们可以从中选取特征波峰或波谷,并将其带入经验公式中。

    在这里插入图片描述
    (3) 如根据上述两幅图,我们分别选择波长为701.5500nm波段与波长为660.8500nm波段。如本文第一部分所示,将其带入文献中所展示的经验比值模型。
    在这里插入图片描述
    在这里插入图片描述

    (4) 选择“Toolbox”→“Model Maker”→“Model Maker”,在弹出的New_Model窗口中配置模型。
    (5) 经过上述ERDAS建模过程,我们得到依据经验比值方法计算的太湖水体叶绿素a含量反演结果图像。可以看到,图像中出现大量负值。这是由于我们直接选用了文献中的参数,且原始图像未经过大气校正造成的。
    在这里插入图片描述
    在这里插入图片描述

    (6) 此外,依据上述近似的思路,利用光谱的一阶微分数值计算太湖水体叶绿素a含量。其中,本次所选择的一阶微分处理特征波段分别为波长为691.3700nm与671.0200nm的两个波段。

    在这里插入图片描述

    (7) 分别将以上两幅经验比值法、一阶微分法计算得出的叶绿素a含量结果制作为专题地图。上述经验比值法计算得到结果存在较多负值,故此处暂不展示其专题地图——大家继续往后看即可。
    在这里插入图片描述

    3 大气校正及经验比值法波段调整

    由以上结果可知,不进行大气校正,所得叶绿素a含量反演结果精度较低,甚至经验比值法计算得到结果存在较多负值,肯定是不对的。因此,这一部分我们基于以下两个方面,对叶绿素a含量反演精度加以提升:

    1. 进行大气校正;
    2. 对出了问题的经验比值法所选用的波段加以调整。

    3.1 转换文件数据格式

    (1) 选择“Basic Tools”→“Convert Data (BSQ,BIL,BIP)”,在弹出的文件选择窗口中选择经过波段合成、编辑头文件操作后的结果图像,点击“OK”。

    在这里插入图片描述
    (2) 调整为“BIL”数据存储格式后,即可使用该图像文件进行FLAASH大气校正。

    3.2 FLAASH大气校正

    (1) 选择“Basic Tools”→“Preprocessing”→“Calibration Utilities”→“FLAASH”,在弹出的转换因子窗口中选择第二项,即单一因子适用于所有波段的情况。由于我们本次所使用数据原有光谱数值为绝对辐射值的标准单位,即(μW)/(cm2nmsr),这一单位为FLAASH方法所能利用的单位,故我们需要将转换因子设定为1.00。

    在这里插入图片描述
    (2) 随后弹出的配置对话框中,首先选择输入图像文件、输出图像文件目录及名称,同时依据遥感影像的元数据,配置其中心点经纬度、传感器类型(传感器类型一旦选定,系统将会自动匹配传感器高度与像元大小这两个参数)、航行时间;同时依据实际研究区的情况,配置平均海拔高度这一选项;其次,选择合适的地球大气模型和气溶胶模型。其中,地球大气模型需要根据一张标准查找表确定,气溶胶模型则依据研究区域实际情况加以确定。设置FLAASH大气校正参数如下。其中,Hyperion产品Level 1数据的航行时间、图像经纬度等元数据可以通过解压缩后的文件夹中“EO1H1190382004095110KZ_SGS_01.fgdc”文件查阅。

    在这里插入图片描述
    (3) 接下来,选择“Multispectral Settings”。当基本设置中设置了水汽反演模型和气溶胶模型时,我们需要相应地将多光谱相关属性加以配置,配置完毕后点击“OK”。
    (4) 全部设置完毕后,可以点击右下角的“Save”按钮,从而将本次对于FLAASH方法的全部属性配置保存;在后期若对同样的遥感数据进行同样的操作时,保存属性配置可以免去多次重复调整相关参数的工作,从而提高效率。
    (5) 点击“Apply”,将会开始执行FLAASH大气校正。但在开始前,将会弹出一个FWHM配置窗口。这是由于,FLAASH若针对高光谱图像加以执行,需要配置这一属性。

    在这里插入图片描述
    (6) 通过与导入中心波长相似的操作,我们将FWHM同样通过Excel软件与TXT文件导入。

    在这里插入图片描述
    (7) 随后开始执行FLAASH大气校正。但发现,FLAASH大气校正执行一定时间后,总会出现未包含错误原因的报错提示,并自动暂停执行过程。

    在这里插入图片描述
    (8) 多次尝试依据网络资源对FLAASH大气校正参数加以调整,但这一问题依然存在。针对这一问题的分析附于下文3.4的(7)部分。

    3.3 QUAC快速大气校正

    (1) FLAASH大气校正无法执行后,尝试使用QUAC快速大气校正方法。选择“Basic Tools”→“Preprocessing”→“Calibration Utilities”→“Quick Atmospheric Correction”,在弹出的文件选择窗口中选择待处理图像文件,点击“OK”。

    在这里插入图片描述
    (2) 得到结果图像文件后,在其上任意位置处右键,选择“Z Profile (Spectrum)”,查看图像中任意位置对应波谱信息。

    在这里插入图片描述

    (3) 将得到的QUAC快速大气校正图像与校正前图像光谱曲线加以对比,可以看到针对地表植被,光谱曲线走势更加符合植被特征曲线,尤其是550nm附近的反射峰、670nm附近的吸收谷、700nm附近的吸收陡坡、900至1300nm附近的平缓趋势等,QUAC大气校正结果图像均有着较好的展现。

    3.4 经验比值法调整

    (1) 通过ENVI软件QUAC快速大气校正后,尝试将大气校正后的结果图像重新带入第一次未成功的经验比值模型中,再一次计算这种方法得到的叶绿素a含量。

    在这里插入图片描述

    在这里插入图片描述
    (2) 大气校正后的结果图像为包括大部分非水体的原有范围区域图像,因此需要对其加以重新裁剪、监督分类。

    在这里插入图片描述
    (3) 本次监督分类时,尝试用不同R、G、B波段配置方式。随后发现,若采用Hyperion产品的真彩色配色方式(29:23:16)得到的结果虽然更加符合实际情况,但其图像像素之间变得较为难以识别,不利于监督分类。
    在这里插入图片描述

    (4) 另一方面,在重新进行监督分类时,发现总是会报出如下错误。多次尝试,均无法避免。
    在这里插入图片描述

    (5) 因此,对输入的QUAC大气校正结果图像加以光谱曲线加以复制,并导出到Excel软件中验证。

    在这里插入图片描述
    (6) 可以看到,上述光谱曲线在波长前一部分并无问题,而在930nm左右之后,其数值均不再发生变化。经过多次重复、尝试,找到问题所在——由于波段合成、编辑头文件、大气校正等过程均在ENVI软件中执行,因此使得输出的文件名越来越长(ENVI软件每完成一步光谱处理均会在原波段名称前增添一个操作名称单词);而过长的ENVI文件名超出了ERDAS软件的文件名识别上限,因此软件将B57波段后的所有波段均识别为同一波段,造成错误。错误的图像如上页最后一幅图;下图则为一幅未经过ENVI大气校正而没有产生这一问题的正常高光谱图像的波段信息。

    在这里插入图片描述
    在这里插入图片描述

    (7) Excel软件中,可以清晰地看出由B79波段开始,所有光谱数值不再发生变化。之所以在B79波段开始出现这一问题,一定是因其为SWIR波段的第一个波段,而SWIR波段相较之V-NIR波段多经过一次ENVI软件的波段合成(因为其额外需要由三组波段合成为一组);多经过一次ENVI操作,其名称会多出一个英文单词,从而出现这一问题。且由此分析,FLAASH之所以无法正常完成,亦是因为我是用ENVI进行图像预处理,使得波段名称过长;而FLAASH运行过程中会生成很多新的临时文件,这些新文件的波段命名方式应亦为在各波段前添加单词,从而导致部分波段名称中相互之间唯一不同的字符被“挤出”软件可识别的名称范围,出现重名,导致FLAASH在运行过程中报错。

    在这里插入图片描述
    (8) 因此,为解决经验比值算法出现的叶绿素a含量为负数的问题,决定只得由第二个角度——特征波段入手。因为我们直接运用文献中推导的经验比值模型参数,且特征波段的选择具有一定主观性;因此特征波段亦会影响最终结果取值。

    在这里插入图片描述
    (9) 多次观察反射率光谱曲线,同时适当结合一阶微分、倒数对数光谱曲线 ,以确定出新的特征波段。

    在这里插入图片描述
    (10) 同样利用ERDAS软件,以期建立模型求出采用新特征波段的叶绿素a含量。最终,确定出701.5500nm与691.3700nm这两个特征波段作为选定的波段。上述特征波段确定的叶绿素a含量数值虽仍有小于0的部分,但已有很大改观。
    在这里插入图片描述
    在这里插入图片描述

    (11) 修正后的经验比值模型求得的太湖水体叶绿素a含量专题地图如图所示(经过指正,以下这幅我的专题地图有误,大家所得结果可能与我的不太一样,依据实际情况判断即可)。同上述第一幅专题地图一致,在制作时需要将图像中太湖周围(即陆地部分)的0值转为NoData。
    在这里插入图片描述在这里插入图片描述

    4 经验比值法、一阶微分法的叶绿素a含量反演结果对比

    (1) 针对上述两种水体叶绿素a反演算法得到的两幅专题地图,对比、分析如下。
    在这里插入图片描述

    (2) 其中,左侧图像为经验比值模型反演结果,右侧图像为一阶微分模型反演结果。
    (3) 在取值范围上,经验比值结果取值分布于0至0.040mg/L左右,一阶微分结果取值分布于0.020至0.065mg/L左右。依据文献及实际情况,可以看到上述两种方法得到的结果整体均偏小,尤其是缺少含量较高的数值(我国东部亚热带、温带内陆湖泊水体叶绿素a含量最大值往往可达0.110至0.120mg/L以上);其中,经验比值得到的结果明显小于一阶微分结果,前者最大值甚至仅仅达到后者最小值左右的水平;且经验比值结果中包含极少量负值像素。由此可知,在数值方面,一阶微分模型精度更高一些。当然,我们本次模型中常量参数均直接利用文献给出的数值,而实际中这些参数的数值是需要依据地表水体实测叶绿素含量、各波段与叶绿素含量的相关系数等加以推导、验证得出的,即上述取值范围有很大一部分是由于参数选取并不是很准确导致的。
    (4) 在分布趋势上,可以看到两种方法得出的结果均符合“岸边含量达到最高,水体中心含量相对最低,由岸边至水体中心含量逐渐降低”这一情况。而不同的是,经验比值得到的结果随离岸距离的变化极其明显,呈现由东北至西南方向的条带状分布;甚至在太湖这一区域的西南部分,由于河岸向水体中心方向明显延申,导致其随亦是近岸区域,但反演得到的叶绿素a含量并没有随之升高——这样的条带状分布过于明显,可能是由于未经处理的光谱曲线含有的噪声或共线性过多,光谱的大气影响等未很好消除而导致的。而一阶微分得到的结果图像相对前者而言,其叶绿素a含量随水体向内延伸而变化的幅度较小,并不是特别明显。而两幅图中,叶绿素a含量最高和最低的区域则较为一致,均为西南圆弧状岸边处含量整体较高,东北部水域含量整体较低。
    (5) 综上所述,可以看出无论是在具体数值取值范围上,还是在数值的空间分布上,一阶微分得到的结果相对较为真实。

    5 其他探究内容

    5.1 基于ENVI实现微分计算

    (1) 在ENVI 5.3 (64-bit) 软件中,可以借助“deriv”函数实现光谱曲线的微分计算。
    在这里插入图片描述

    (2) 与经典版不同的是,ENVI 5.3 (64-bit) 软件在编辑头文件时有所变化,但整体步骤流程一致。此外,太湖区域的提取亦可以在ENVI软件中通过剪裁实现。
    在这里插入图片描述
    在这里插入图片描述

    5.2 水体表层沉积物平均粒径反演

    (1) 依据本文第一部分原理,反演计算太湖区域水体表层沉积物平均粒径。

    在这里插入图片描述

    (2) 求出专题地图。

    在这里插入图片描述
    (3) 由专题地图可以看出,所求得的水体表层沉积物平均粒径取值主要集中于6.0Ф以上;其结果图像中高值分布区域与前期在做监督分类时采用的假彩色图像中太湖水体浑浊(即泛白色)位置(如下图所示)十分一致,可见原本遥感图像中浑浊其实为水体中表层沉积物;沉积物平均粒径越大,水体越浑浊。而上述结果中的条纹状现象可能为特征波段选取并不是很恰当,其间具有部分共线性导致的。
    在这里插入图片描述

    5.3 土壤有机质反演

    (1) 依据文献,利用两个特征波段各自倒数对数的一阶微分的比值与土壤有机质含量的二次回归方程,经过同样的建模、制图过程,得出太湖周边地区土壤有机质反演专题地图。
    (2) 从中可以看出,土壤有机质含量较高的区域主要集中于城镇中心、道路两旁,这可能是由于人为、交通等导致的碳积累;大部分地区有机质含量处于0.5%至1.0%左右,这个值或许有些小于常见水平。图像中的小面积空白区域多为非太湖湖泊、水田、鱼塘等水体。

    在这里插入图片描述

    在这里插入图片描述

    本文公式、反演算法等可以参考的文献

    在这里插入图片描述

    欢迎关注公众号:疯狂学习GIS
    在这里插入图片描述

    展开全文
  • 基于二次时域微分解析的油纸绝缘介质响应参数辨识
  • 在具有明确的作用方程式的情况下,直接用微分法可以进行误差分析和计算。例如图1所示的自准直仪测角误差中,除原理误差外,还存在制造误差。那么,自准直仪的测角误差与哪些零件的制造误差有关呢?由式(4-2)得 ...
  • §8.6 微分法在几何上的应用 一、空间曲线的切线与法平面 1、曲线由参数方程给出的情形 设空间曲线的参数方程为  (1) 假定(1)式中的三个函数均可导。 考虑上对应于的一点及对应于的邻近一点,其割线的方程...

    §8.6  微分法在几何上的应用

    一、空间曲线的切线与法平面

    1、曲线由参数方程给出的情形

    设空间曲线的参数方程为

                                 (1)

    假定(1)式中的三个函数均可导。

    考虑上对应于的一点及对应于的邻近一点,其割线的方程为

    对等式同除以

    时,,曲线在点处的切线方程为

                                  (2)

    这里自然假定了  不能都为零。

    切线的方向向量称为曲线的切向量,向量

    就是曲线在点处的一个切向量。

    过点与切线垂直的平面称为曲线在点处的法平面,它是过点,以为法向量的平面,此法平面方程为

             (3)

    2、曲线由特殊参数方程给出的情形

      此方程可看作 

    处可导,则,曲线在点处的切线方程为

                                      (4)

    曲线在点处的法平面方程为

                  (5)

    3、曲线由一般方程给出的情形

    是曲线上的一点,此函数方程组可确定的隐函数,即曲线可用(隐式)方程  来表示。

    由第2部分的讨论,现在的关键是求

    看作的隐函数,方程两边分别对求导数,可得

         ð   

     

    ð 

    ð   

    ð     ð 

    类似地,有

    曲线在点处的切向量本来为,但也可取向量

    即 

    曲线的切线方程为

                             (6)

    曲线的法平面方程为

       (7)

    当然,上述推导需要一些条件,具有一阶连续偏导数,且

    中至少有一个不为零。

    【例1】求曲线

    在点处的切线方程与法平面方程。

    解:

     

     

    曲线的切线方程为  

    曲线的法平面方程为 

    二、曲面的切平面与法线

    1、曲面方程由给出的情形

    设曲面由方程

                                             (9)

    给出,上的一点,假设函数偏导数在该点连续且不同时为零

    上,过点任意引一条曲线,设它的参数方程为

    对应于参数,且不全为零。

    则曲线在点的切线方程为

    下证事实:

    上过点且具有切线的任何曲线 ,它们在点处的切线均位于同一平面。

    因为曲线在曲面上,故有

    据假设有  ,即

         (10)

    引入向量

    (10)式表明:

    因为是过点且在上的任意一条曲线,它们在点的切线均垂直于同一非零向量,所以,上过点的一切曲线在点的切线都位于同一个平面上。

    这个平面称为曲面在点切平面,其切平面方程为

           (11)

    过点而垂直于切平面(11)的直线称为曲面在该点的法线,其法线方程为:

                    (12)

    曲面在一点的切平面之法向量称为曲面在该点的法向量,因此,向量

    便是曲面在点处的一个法向量。

    2、曲面方程由给出的情形

    若曲面由方程

    给出,令

    则 

    当偏导数在点连续时,曲面在点的切平面方程为

                                 (14)

    曲面的法向量有两个

    对于第一式,法向量的方向余弦为

    ,法向量与轴正向的夹角应为锐角,故此法向量的指向是朝上的。自然地,另一个法向量的指向是朝下的

    (13)式具有鲜明的几何意义

    方程的右端恰好是函数在点处的全微分;

    方程的左端是切平面上点的竖坐标的增量。

    特别地,当  时

    曲面在点处的切平面为,此切平面平行于坐标面,即曲面在点处具有水平的切平面

    【例2】求球面在点处的切平面及法线方程。

    解:

    切平面方程为

    法线方程为

    因为点在法线上,可见法线通过球心。

     

     

    展开全文
  • 微分零交叉是目前应用最多的算法之一, 但是在信噪比较低情况下会造成较多误判。根据激光雷达回波信号的特点, 对微分零交叉进行了改进。改进的算法中参考了候选点信号强度特性和前后时刻的云层信息, 有效地修正了...
  • 在具有明确的作用方程式的情况下,直接用微分法可以进行误差分析和计算。例如图1所示的自准直仪测角误差中,除原理误差外,还存在制造误差。那么,自准直仪的测角误差与哪些零件的制造误差有关呢?由式(4-2)得 ...
  • 圆的两种生成算法(角度微分法、Bresenham算法) 文章目录1.角度微分法的原理2.角度微分法的实现(基于matlab)3.Bresenham 算法的原理4.Bresenham 算法的实现(基于matlab) 1.角度微分法的原理 圆的角度微分法是...

    圆的两种生成算法(角度微分法、Bresenham算法)





    1.角度微分法的原理

    圆的角度微分法是用圆的内接正多边形来逼近该圆。
    若我们设圆的参数方程为:
    Alt
    其中𝜑为旋转角,0° ≤ φ ≤ 360° 现把该圆 n 等分,用 n 条直线来逐次逼近该圆,则各直线的端点坐标计算如下:
    设旋转角𝜑的起始角、终止角分别为α、β,且满足关系式0° ≤ α ≤ β ≤ 360°
    把圆进行 n 等分,其对应的角度增量为:𝛥𝜑 = 2𝜋 ∕ 𝑛
    则可以得到第 i 点所对应的角度与端点坐标为(其中 i = 0 , 1 , 2 , …… , n):
    Alt
    上式即为生成圆的基本公式之一。我们不难看出,n 越大,精度越高,但计算也越 复杂。下一步是如何确定 R 和 n 的关系?
    可以这样认为:如果圆的内接多边形的边线与理想圆弧之间的距离在一个刻度单位之内时,其近似显示效果是能接受的。如下图所示:
    Alt
    即:如果线段 ab 的距离在一个刻度单位内都是能被接受的,即要求 ab < 1
    那 ab 的距离等于多少呢?显然 𝑎𝑏 = 𝑅 − 𝑅 ⋅ 𝑐𝑜𝑠𝛥𝜑 < 1
    对于𝑐𝑜𝑠𝛥𝜑,可按泰勒级数展开取前两项,得到𝑐𝑜𝑠𝛥𝜑 = 1 − 𝛥𝜑2 / 2
    将其带入 ab 中即得到 𝑅 − 𝑅 ⋅ (1 − 𝛥𝜑2 / 2 ) < 1,解之可得𝛥𝜑 < √(2 / 𝑅)
    将𝛥𝜑 = 2𝜋 / 𝑛带入上式可得𝑛 > 𝜋 ⋅ √(2𝑅)
    用经验公式取代得到 𝑛 = 0.1𝑅 + 15
    这就是圆的半径与其内接正多边形边数之间的关系




    2.角度微分法的实现(基于matlab)

    下面给出利用matlab实现该算法的完整代码:

    function main         						%测试主函数
    clc
    clear;
    DifferentiationOfCircle(100,'r')			%调用角度微分法
    
    function DifferentiationOfCircle(R,color)   %圆的角度微分算法
        n = 0.1*R+15;                           %计算出需要微分的次数
        delta = 2*pi/n;                         %从而得到每次需要增加的角度
        alpha = delta;                          %初始化alpha
        x0 = R,y0 = 0;                          %初始化x0,y0
        hold on
        while(n>0)
            x = R*cos(alpha);
            y = R*sin(alpha);
            plot([x0,x],[y0,y]);                %绘制点(x0,y0)(x,y)之间的线段
            alpha = alpha + delta;              %更新alpha
            x0 = x , y0 = y;                    %更新x0、y0
            n = n-1;                            %n减1
        end
        grid minor
        hold off
    

    在matlab中运行上述代码,实验结果如下:
    Alt
    可以看出,角度微分法在绘圆时具有明显的分段现象。当然,我们可以通过加大 n 的值来适当减轻分段效果。不过这样的代价是会使得程序需花费更多的时间来进行微分。 所以这才有了𝑛 = 0.1𝑅 + 15 这一公式的出现。




    3.Bresenham 算法的原理

    圆可被定义为到给定中心位置 O(圆心)距离为 r 的点集。圆心位于原点的圆有四条对称轴 x=0,y=0, x=y 和 x=-y。若已知圆弧上一点(x,y),可以得到其关于四条对称轴的其它 7 个点,这种性质称为八分对称性(如下图所示)。因此,只要扫描转换八分之一圆弧,就可以求出整个圆弧的象素集。
    Alt
    所以接下来我们的重心就放在:如何将某段圆弧尽可能描绘清楚
    对于上图中,点(y,x)所在那一段圆弧而言,其特征是 x 坐标变化率大于 y 坐标变化 率,且随着 x 的增加,y 在减小。因此逼近该圆弧的绘点步进规律是:每推算一次,其 x 坐标都加 1,而 y 坐标根据相应误差决定是否减 1。我们易知该圆中最上方的那个点的 坐标为(0,R),将其设为(𝑥0,𝑦0),则有ⅆ(𝑥0,𝑦0) = 0)。对于下一点而言,有两种可能:
    ① (𝑥0+1,𝑦0)
    ② (𝑥0+1,𝑦0−1)
    现在我们的任务就是确定到底取谁。
    在高中时,我们对于圆的定义式为:ⅆ(𝑥,𝑦) = 𝑥2 + 𝑦2 − 𝑅2
    这样定义的意义在于:对于某个给定点(𝑥0,𝑦0),其值的正负能反应该点位于圆内 还是圆外。具体的规则如下:
    若ⅆ(𝑥0,𝑦0) < 0,则表示该点位于圆内;
    若ⅆ(𝑥0,𝑦0) > 0,则表示该点位于圆外;
    若ⅆ(𝑥0,𝑦0) = 0,则表示该点位于圆内;
    在确定到底取①还是②时,我们实际上并不关心它们到底在圆内还是圆外,我们只 关心它们到底谁距离圆的那一段圆弧更近。因此,我们可以将这两个点带入圆的定义式 中,即得到ⅆ(𝑥0+1,𝑦0)、ⅆ(𝑥0+1,𝑦0−1),通过加绝对值的方式来获得这两点距离圆弧的绝对距离,最后取其中的较小值作为下一个点的位置即可。

    因此,如果我们设:
    ⅆ(𝑈𝑖) = (𝑥𝑖−1 + 1)2 + (𝑦𝑖)2 − 𝑅2
    ⅆ(𝐷𝑖) = (𝑥𝑖−1 + 1)2 + (𝑦𝑖 − 1)2 − 𝑅2
    则可以将圆的 Bresenham 算法的偏差判别式定义为:

    ⅆ𝑖 = |ⅆ(𝑈𝑖)| − |ⅆ(𝐷𝑖)|

    其原则是:
    ① 当 ⅆ𝑖 < 0 时,选 𝑈𝑖 作为下一个绘制的像素点
    ② 当 ⅆ𝑖 ≥ 0 时,选 𝐷𝑖 作为下一个绘制的像素点
    将上面偏差判别式的绝对值去掉,得到:

    ⅆ𝑖 = ⅆ(𝑈𝑖) + ⅆ(𝐷𝑖)

    根据这个判别式,我们可通过一个循环来将八分之一圆弧绘制出,同时通过坐标变 换以将整个圆画出。但是,这个式子里的运算偏多,我们需要对其进行优化,以使整个 绘制过程中,尽量是以迭代运算为主,从而减少复杂的乘幂运算次数。

    首先是将𝑈𝑖(𝑥𝑖−1 + 1, 𝑦𝑖−1),𝐷𝑖(𝑥𝑖−1 + 1,𝑦𝑖−1 − 1)的坐标带入得到:
    𝑖 = 2(𝑥𝑖−1 + 1)2 + (𝑦𝑖−1 − 1)2 + (𝑦𝑖−1)2 − 2𝑅2

    然后增加下标值可得到:
    𝑖+1 = 2(𝑥𝑖 + 1)2 + (𝑦𝑖 − 1)2 + 𝑦𝑖2 − 2𝑅2

    将其作差(ⅆ𝑖+1 − ⅆ𝑖),得到:
    𝛥ⅆ = 2((𝑥𝑖 + 1)2 − (𝑥𝑖−1 + 1)2) + ((𝑦𝑖 − 1)2 − (𝑦𝑖−1 − 1)2) + (𝑦𝑖2 − (𝑦𝑖−1)2)

    ① 当ⅆ𝑖 < 0时,选择的下一点是𝑈𝑖,此时𝑥𝑖 = 𝑥𝑖−1 + 1, 𝑦𝑖 = 𝑦𝑖−1,带入上式可得: 𝛥ⅆ = 4𝑥𝑖−1 + 6;
    ② 当ⅆ𝑖 ≥ 0时,选择的下一点是𝐷𝑖,此时𝑥𝑖 = 𝑥𝑖−1 + 1,𝑦𝑖 = 𝑦𝑖−1 − 1,带入上式 可得: 𝛥ⅆ = 4(𝑥𝑖−1 − 𝑦𝑖−1) + 10 。
    这样一来,对于 ⅆ𝑖 的求解就由之前的算术求值变为了迭代求值,大大减少了程序的 运算时间。




    4.Bresenham 算法的实现(基于matlab)

    下面给出利用matlab实现该算法的完整代码:

    function main         	%测试函数
    clc
    clear;
    Bresenham(100)			%调用Bresenham算法
    
    function Plot(x,y)      %此函数会将某个点集在对应的8个区域中所构成的弧线绘制出,从而构成一个圆
        plot(x,y,'r')       %用红色的线绘制第一象限的圆弧
        plot(y,x,'r')
        plot(-x,y,'b')      %用蓝色的线绘制第二象限的圆弧
        plot(-y,x,'b')
        plot(-x,-y','y')    %用黄色的线绘制第三象限的圆弧
        plot(-y,-x,'y')
        plot(x,-y,'g')      %用绿色的线绘制第四象限的圆弧
        plot(y,-x','g')
        
    function Bresenham(R)         %Bresenham算法
        d = 0 , n = 1;            %初始情况下d=0
        x(1) = 0 , y(1) = R;
        while(x < y)
            if(d < 0)
                d = d + 4*x(n) + 6;
                y(n+1) = y(n);
            else
                d = d + 4*(x(n)-y(n)) + 10;
                y(n+1) = y(n) - 1;
            end
            x(n+1) = x(n) + 1;
            n = n + 1;
        end
        hold on;
        Plot(x,y);
        grid minor;
        hold off;
    

    在matlab中运行上述代码,实验结果如下:
    Alt
    可以看出,Bresenham 算法在绘圆时不会出现明显的分段现象,但是却存在“抖动”弧。出现这一现象的原因是由于 plot 函数是对像素点进行连结绘制,这个函数并不具有平滑拟合的效果。但如果仅对此算法在描绘圆弧的点上进行评估的话,该算法的效果明显是优于微分算法的。




    展开全文
  • 第九章、多元函数微分法及其应用 知识逻辑结构图 考研考试内容 多元函数的概念,二元函数的几何意义,二元函数的极限(极限存在的判定)和连续的概念,有界闭区域上多元连续函数的性质(有界性,最值存在,介值定理...
  • 这篇文章翻译自官方教程Automatic Derivatives并且参考了少年的此间的博客文章...在这一章中,我们将考虑几种不同的方法来在这些特殊情况下使用自动微分法。 现在我们考虑一个优化问题。寻找参数θθ\theta和ttt...
  • 第八章 多元函数微分法及其应用

    千次阅读 2014-03-25 14:49:49
    第七讲 方向导数与梯度 教学目的 使学生理解方向导数与梯度的概念,掌握方向导数与...是与同方向的单位向量射线的参数方程为   设函数在点的某个邻域内有定义,为上另一点,且,到的距离  若 当沿着趋于即
  • 论文研究-微分建模及其在疲劳可靠性研究中的应用.pdf, 很多工程问题可用微分方程描述,由此导出了微分建模方法.该方法的实质是直接用微分方程来描述客观物理现象,同时...
  • 这篇文章翻译自官方教程Analytic Derivatives并且参考了少年的此间...待确定参数方程如下: y=b1(1+eb2−b3x)1/b4y=b1(1+eb2−b3x)1/b4y = \frac{b_1}{(1+e^{b_2-b_3x})^{1/b_4}} 现在给定一系列的对应数据点{xi,y...
  • 一阶线性微分方程的初等积分变量分离方程与变量变换可分离变量的微分方程齐次方程可化为齐次方程的类型线性方程与常数变易法一阶线性齐次微分方程一阶线性非齐次微分方程伯努利方程黎卡提方程恰当方程与积分因子...
  • 一阶线性微分方程的初等积分变量分离方程与变量变换可分离变量的微分方程例子齐次方程例子可化为齐次方程的类型例子线性方程与常数变易法一阶线性齐次微分方程一阶线性非齐次微分方程伯努利方程黎卡提方程恰当方程...
  • 自动微分

    2019-12-10 17:32:36
      梯度下降用于求损失函数的最优值,梯度下降是通过计算参数与损失函数的梯度并在梯度的方向不断迭代求得极值;但是在机器学习、深度学习中很多求导往往是很复杂的,手动使用链式法则求导很容易出错,借助于...
  • 平均的两个步骤: 一是假设解的形式,这个形式由非线性激励和/或非线性阻尼中的小参数取零时确定,即系统在无阻尼无激励时的精确解形式,然后令幅值和相位变为时间的函数 二是认为幅值和相位的变化是小量,在一...
  • Lane-Emden型方程是非线性微分方程,出现在许多领域,例如恒星结构,放射性冷却和星系团建模。 在这项工作中,该方程是使用一种称为辅助参数的半参数方法的半分析方法来研究的。 在应用技术中,将未知的辅助参数插入...
  • matlab 龙格-库塔 求解常微分方程

    万次阅读 多人点赞 2012-06-18 11:29:19
    最近学习分室模型,里面碰到了用matlab 龙格...function [x,y]=runge_kutta1(ufunc,y0,h,a,b)%参数表顺序依次是微分方程组的函数名称,初始值向量,步长,时间起点,时间终点(参数形式参考了ode45函数) n=floor((b-a
  • 前在经典控制理论中广泛使用的分析设计方法--频率和根轨迹,不是直接求解微分方程,而是采用与微分方程有关的另一种数学模型——传递函数,间接地分析系统结构参数对响应的影响。所以传递函数是一个极其重要的...
  • 针对一类由半线性抛物型偏微分方程描述的分布参数系统,考虑系统边界为Robin或混合边界条件的情况,提出基于边界控制的控制策略,并研究其镇定问题.首先,根据无穷维抽象发展方程理论和Lumer-Phillips理论证明闭环系统的...
  • 如果我们将x理解成delta_x*n,则u(x,t)可以看成只是t的函数,n作为参数来控制x,即定义: 这样我们就可以把这个偏微分方程看做是一个常微分方程 对于常微分方程的C++实现,可以先看看...
  • 参数函数求导 1. 利用莱布尼茨定理求高阶导 只看两点: 1、常用导数的高阶公式 2、例题 例题: 2.隐函数求导 这种方程里面y是x的函数,但是不显性。 例题1 设 y=y(x),y2−2xy+9=0y = y(x), y^2 - 2xy +9 = 0y=y(x)...
  • 经典信号建模(classical modeling method for signal)前面已经指出,医学信号处理的目的是提取包含于随机信号中的确定性成分,以便在一定的准确性(最小二乘意义)上进行预测。这就是建立各种各样的确定性数学...
  • 关于一元函数微分学,专升本数学考试要求包括:(一)导数与微分1.理解导数的概念及其几何意义,了解左导数与右导数的定义,理解函数的可导性与连续性的关系,会用定义求函数在...掌握对数求导参数方程求导。...
  • 这里对使用python求解常微分方程提供两种思路,一种是自己编程实现欧拉,改进欧拉或者四阶龙格库塔,这样有助于理解上述三种数值计算方法的原理;一种是调用python已有的库,不再重复造轮子。本文对上述两种思路...
  • 01微分几何

    2020-07-11 23:27:58
    微分几何 曲线(Curve) 曲线是一个函数的图像表示。就比如高中物理中最经常出现的位移-时间函数,他在二维平面由s-t表示,也就变成了参数方程。曲线的正切向量定义为曲线的求导,也就是曲线的斜率。把这个向量旋转90...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 429
精华内容 171
关键字:

参数微分法