CN102032903B - 基于Landsat数据源的珊瑚岛礁遥感信息自动提取方法 - Google Patents

基于Landsat数据源的珊瑚岛礁遥感信息自动提取方法 Download PDF

Info

Publication number
CN102032903B
CN102032903B CN2010105675077A CN201010567507A CN102032903B CN 102032903 B CN102032903 B CN 102032903B CN 2010105675077 A CN2010105675077 A CN 2010105675077A CN 201010567507 A CN201010567507 A CN 201010567507A CN 102032903 B CN102032903 B CN 102032903B
Authority
CN
China
Prior art keywords
idx
highland
remote sensing
coral
image
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN2010105675077A
Other languages
English (en)
Other versions
CN102032903A (zh
Inventor
李满春
刘永学
程亮
黄秋昊
江冲亚
张昱
刘成明
李�真
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Nanjing University
Original Assignee
Nanjing University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Nanjing University filed Critical Nanjing University
Priority to CN2010105675077A priority Critical patent/CN102032903B/zh
Publication of CN102032903A publication Critical patent/CN102032903A/zh
Application granted granted Critical
Publication of CN102032903B publication Critical patent/CN102032903B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Image Processing (AREA)

Abstract

本发明公开了基于Landsat数据源的珊瑚岛礁遥感信息自动提取方法，属于遥感影像地物自动提取领域。其步骤为：遥感影像辐射定标；基于图像的遥感影像大气校正；对遥感影像进行MeanShift滤波，消除条带和椒盐噪声；对影像第5波段进行直方图阈值分割，区分非高潮高地和高潮高地；对于高潮高地区分陆地植被和灰沙/建筑；对于非高潮高地，区分低潮高地和非干出质底；对于低潮高地区分海藻/海草和浅水珊瑚；对于非干出质底区分深水珊瑚和开放水体。本发明能够分阶段逐层次分解地提取出陆地植被、灰沙/建筑、海藻/海草、浅水珊瑚、深水珊瑚5种珊瑚岛礁覆盖类型及开放水体背景信息，实验表明，本发明提取精度较高，能够有效改善珊瑚岛礁遥感调查与监测的效率。

Description

基于Landsat数据源的珊瑚岛礁遥感信息自动提取方法
技术领域
本发明涉及一种遥感影像珊瑚岛礁自动提取方法，特别是涉及一种采用分阶段逐层次分解的方式建立了珊瑚岛礁遥感信息提取方法。
背景技术
珊瑚礁是生物多样性最高、最具生态价值的海洋生态系统，具有极高的经济价值、药用价值、旅游价值和科研价值。全面、快速、动态地开展珊瑚岛礁调查是一项非常重要的基础工作。传统的实地测量方法不仅费时费力，而且无法获得大范围的观测数据，特别是我国南海诸岛由于地理位置及国际形势的原因，开展实地调查尤为困难。遥感技术作为大范围的、实时的全球观测技术，是调查和监测大范围珊瑚礁状况的一种有效手段，通过运用不同类型的传感器，能够完成珊瑚礁群落分类、珊瑚礁群落变化监测、珊瑚礁白化监测、珊瑚礁水深测量等任务。目前，珊瑚岛礁遥感信息提取主要是借助目视解译方法和监督分类两种方式，而要满足大范围珊瑚岛礁动态监测的需求，遥感信息提取方法的自动化程度尚有待提高。
潘春梅等人于2002年在期刊《国土资源遥感》总第52期中发表“南沙群岛岛礁地形地貌TM影像特征分析”一文，具体根据南沙海区的地形地貌特征，利用Landsat卫星TM遥感图像，对不同的珊瑚礁、珊瑚沙发育状况和岛礁的地质特征与生态特征进行了分析，研究了水下礁、滩的物质组份、发育程度及TM影像特征，选择TM4、TM2、TM1波段图像组合，提出了岛礁分类方法，然而该方法仅仅提出理论参考和经验方法，难以正在实现全自动化岛礁提取处理。王小龙等于2005年在期刊《海洋科学进展》23卷第4期发表“基于光学遥感的海岛潮间带和湿地信息提取——以东沙岛(礁)为例”论文，利用多源卫星遥感数据和东沙岛水深地形图，结合东沙岛潮汐调和常数，确定潮间带和湿地的范围，用于提取潮间带和湿地的面积，该方法从一定角度上解决了岛礁边界划分以及分类的问题，然而，该方法仍然是基于人工的方法，无法实现大范围大数据量下的岛礁自动提取。栗敏光等于2010年在期刊《国土资源遥感》总第83期中发表“多源多时相遥感数据面向对象海岛识别方法探讨”论文，针对目视解译方法的不足，采用面向对象的自动分类方法实施海岛专题信息提取，并对多源多级遥感图像数据和多时相高分辨率光学图像数据，分别采取了不同的小目标识别与伪信息剔除策略，然而，论文仅仅提出一个探索性方法，以较小的实验区进行初步的实验，还无法真正应用于大范围下的珊瑚岛礁调查。
总体而言，珊瑚岛礁遥感调查还处于探索阶段，以利用TM、ETM+、SPOT等卫星遥感影像的手工、目视解译珊瑚岛礁为主。遥感影像处理过程中，以常规的商业化通用软件为主要工具手段，而具体珊瑚岛礁解译的算法精度、速度、质量及自动化程度还缺少相应的研究，研究工作还处于较低水平的重复劳动阶段，珊瑚岛礁遥感识别有待改善。
发明内容
1.发明要解决的技术问题
本发明的目的是提供基于Landsat数据源的珊瑚岛礁遥感信息自动提取方法，快速有效地提取出陆地植被、灰沙/建筑、海藻/海草、浅水珊瑚、深水珊瑚5种珊瑚岛礁覆盖类型及开放水体背景信息。
2.本发明的技术方案如下：
原理：首先对遥感影像进行正规化处理，减小不同时相、不同轨道的遥感影像数据的光谱特征差异，正规化处理包括辐射定标和基于图像的遥感影像大气校正。然后对正规化处理后的影像进行Mean Shift空间滤波，以减轻条带和椒盐噪声。对影像第5波段进行直方图阈值分割，区分非高潮高地和高潮高地；对于高潮高地，利用NDVI进行直方图阈值分割，区分陆地植被和灰沙/建筑；对于非高潮高地，利用影像第3波段进行直方图阈值分割，区分低潮高地和非干出质底；对于低潮高地，利用NDVI进行直方图阈值分割，区分海藻/海草和浅水珊瑚；对于非干出质底，利用影像第1波段进行直方图阈值分割，区分深水珊瑚和开放水体。直方图阈值分割采用OTSU方法。
基于Landsat数据源的珊瑚岛礁遥感信息自动提取方法，包括以下步骤：
步骤1：对原始遥感影像进行辐射定标，将DN值转为辐射亮度值，定标公式为L＝DN×gain+bias，增益系数gain和偏移量bias采用美国地质调查局(United States Geological Survey，简称USGS)提供的同期参数；
步骤2：对定标后的遥感影像进行大气校正，将辐射亮度值转化为地表反射率，校正公式为：
R = ( L - L p ) · π ( E 0 · cos θ z · T z + E d ) · T v
其中，R为地表反射率，L为步骤1得出的辐射亮度值，E0为大气顶端的太阳辐照度，可由通过太阳常数和日地距离求得，即E0＝Es/D2，θz为太阳天顶角，LP为大气的程辐射能量值，Tz为从太阳到地表入射方向上的大气透射率，Tv为从地表到传感器反射方向上的大气透射率，Ed为经散射作用到达地表的太阳辐照度，具体参数的计算公式为：
Tz＝cosθz
Tv＝cosθv
Ed＝πLp
Lp＝L-0.02·(E0·cosθz·Tz+Ed)·Tv
其中θv为传感器的天顶角；
步骤3：对正规化处理后的影像进行Mean Shift空间滤波，以减轻条带和椒盐噪声；
步骤4：对影像第5波段进行OTSU直方图阈值分割，区分非高潮高地和高潮高地；
步骤5：对于高潮高地，利用NDVI进行OTSU直方图阈值分割，区分陆地植被和灰沙/建筑；
步骤6：对于非高潮高地，利用影像第3波段进行直方图阈值分割，区分低潮高地和非干出质底；
步骤7：对于低潮高地，利用NDVI进行直方图阈值分割，区分海藻/海草和浅水珊瑚；
步骤8：对于非干出质底，利用影像第1波段进行直方图阈值分割，区分深水珊瑚和开放水体。
3.有益效果
与现有技术相比，本发明实现了基于Landsat TM/ETM+数据源的珊瑚岛礁遥感信息自动提取方法，快速有效地提取出陆地植被、灰沙/建筑、海藻/海草、浅水珊瑚、深水珊瑚5种珊瑚岛礁覆盖类型及开放水体背景信息，与目视解译和监督分类方法相比自动化程度大为提高。具体有益效果如下：
第一，本发明针对经正规化处理后的TM/ETM+遥感影像的特点，能够分阶段逐层次分解地提取出陆地植被、灰沙/建筑、海藻/海草、浅水珊瑚、深水珊瑚5种珊瑚岛礁覆盖类型及开放水体背景信息，整个过程无需输入任何参数或人工干预。
第二，本发明提取精度较高，能够有效改善珊瑚岛礁遥感调查与监测的效率。
附图说明
图1为珊瑚岛礁遥感信息提取模型；
图2为试验区遥感影像(431波段组合，经拉伸处理)，其中(a)东沙岛；(b)永兴岛；(c)杨信沙洲；(d)司令礁；(e)永暑礁(西南部)；(f)美济礁；
图3为试验区珊瑚岛礁遥感信息自动提取结果，其中(a)东沙岛；(b)永兴岛；(c)杨信沙洲；(d)司令礁；(e)永暑礁(西南部)；(f)美济礁。
具体实施方式
采用附图2所示遥感影像作为待提取珊瑚岛礁影像。采用C++编程语言在Visual Studio 2005平台下实现本方法的6个步骤，遥感影像数据的读写操作通过开源地理数据格式转换类库GDAL 1.60实现。
步骤1：以GDAL为影像数据读写工具，利用利用GDAL.Open方法读取影像灰度值至大小与原始遥感影像一致的整型数组，遍历数组，按照公式L＝DN×gain+bias将影像DN值转化为辐射亮度值。增益系数gain和偏移量bias采用USGS提供的同期参数；
步骤2：对步骤1中定标后的遥感影像进行大气校正，将辐射亮度值转化为地表反射率，校正公式为：
R = ( L - L p ) · π ( E 0 · cos θ z · T z + E d ) · T v
其中，R为地表反射率，L为步骤1得出的辐射亮度值，E0为大气顶端的太阳辐照度，可由通过太阳常数和日地距离求得，即E0＝Es/D2，θz为太阳天顶角，LP为大气的程辐射能量值，Tz为从太阳到地表入射方向上的大气透射率，Tv为从地表到传感器反射方向上的大气透射率，Ed为经散射作用到达地表的太阳辐照度，具体参数Tz、Tv、Ed、LP的计算公式为：
Tz＝cosθz
Tv＝cosθv
Ed＝πLp
Lp＝L-0.02·(E0·cosθz·Tz+Ed)·Tv
其中，θv为传感器的天顶角；
步骤3：对经步骤2大气校正后的影像进行Mean Shift空间滤波，以减轻条带和椒盐噪声，得到平滑后的TM反射率图像。实现方法的核心代码为：
    void SpatialclusterBasedSegmentAlgorithm::meanShiftSegment(unsigned
char*src，unsigned char*des，int w，int h，int dms，int band1，int band2，int band3，
double winRadius，double dnRadius，int nfilterIter，double dnMinDis，int nclusters，int
nclusterIter)
    {
        //定义变量
        theInputData＝new double[m_size*m_dms]；
        theOutputData＝new int[m_size*m_dms]；
        for(idx＝0；idx＜m_size*m_dms；idx++)
        {
            theInputData[idx]＝src[idx]；
        }
        //滤波
        thefilteredData＝new double[m_size*m_dms]；
        filter(theInputData，thefilteredData，m_band1，m_band2，m_band3，m_winRadiu
s，m_dnRadius，m_nfilterIter)；
        //分割
        thesegmentedData＝new double[m_size*m_dms]；
        labels＝new long[m_size]；
    segment(thefilteredData，thesegmentedData，&cntObjectData，m_dms，labels，m
_dnMinDis)；
        //拓扑
        adjacents＝new long*[cntObjectData]；
        cntMaxAdacents＝(long)sqrt((double)cntObjectData)；
        for(idxObject＝0；idxObject＜cntObjectData；idxObject++)
        {
            adjacents[idxObject]＝new long[cntMaxAdacents]；
        }
        cntAdacents＝new long[cntObjectData]；
        buildTopology(labels，cntObjectData，adjacents，cntAdacents)；
        //聚类
        theclusteredData＝new double[cntObjectData]；
        cluster(thesegmentedData，theclusteredData，cntObjectData，m_dms，adjacents，
cntAdacents，m_nclusters，m_nclusterIter)；
        //对象-＞像元
        theClassData＝new double[m_size]；
        for(idx＝0；idx＜m_size；idx++)
        {
            idxObject＝labels[idx]；
            theClassData[idx]＝theclusteredData[idxObject]；
        }
        //分割
        thesegmentedData＝new double[m_size*m_dms]；
        labels＝new long[m_size]；
        segment(theClassData，thesegmentedData，&cntObj ectData，1，labels，0)；
        //赋值
        long*cntPixels＝new long[cntObjectData]；
        for(idxObject＝0；idxObject＜cntObjectData；idxObject++)
        {
            cntPixels[idxObject]＝0；
        }
        for(idx＝0；idx＜m_size；idx++)
        {
            //当前像元所属对象的编号
            idxObject＝labels[idx]；
            cntPixels[idxObject]++；
            //遍历特征，累加
            for(k＝0；k＜m_dms；k++)
                thesegmentedData[idxObject+k*cntObjectData]+＝
theInputData[idx+k*m_size]；
        }
        //均值
        for(idxObject＝0；idxObject＜cntObjectData；idxObject++)
        {
            for(k＝0；k＜m_dms；k++)
                thesegmentedData[idxObject+k*cntObjectData]/＝
cntPixels[idxObject]；
        }
        for(idx＝0；idx＜m_size；idx++)
        {
            idxObject＝labels[idx]；
            for(k＝0；k＜m_dms；k++)
                theOutputData[idx+k*m_size]＝
(int)thesegmentedData[idxObject+k*cntObjectData]；
        }
        for(idx＝0；idx＜m_size*m_dms；idx++)
            des[idx]＝theOutputData[idx]；
    }
    void SpatialclusterBasedSegmentAlgorithm::filter(double*theInputData，
double*theOutputData，int band1，int band2，int band3，double hs，double hr，int
maxIter)
    {
        int i，j，k；
        long idx；
        double*sdata；
        long idxs；
        int idxBucFt；
        double sMins；
        double sMaxs；
        long idxd；
        //路径点特征向量
        double*yk＝new double[m_dms+2]；
        //均值偏移特征向量
        double*Mh＝new double[m_dms+2]；
double wsum；
double dis；
bool flg；
int nIter；
double w；
double disShift；
int dimension；
int ii，jj，p，q；
long idxNei；
int component1，component2，component3；
component1＝band1-1；
component2＝band2-1；
component3＝band3-1；
double*ykNei；
ykNei＝new double[m_dms+2]；
double*shift；
shift＝new double[m_dms+2]；
double diff1，diff2，diff3；
//遍历所有像元
for(idx＝0；idx＜m_size；idx++)
{
    i＝idx/m_width；
    j＝idx％m_width；
    //特征向量
    yk[0]＝i；
    yk[1]＝j；
    for(k＝0；k＜m_dms；k++)
        yk[k+2]＝theInputData[idx+k*m_size]；
nIter＝0；
flg＝true；
while(flg)
{
    for(k＝0；k＜m_dms+2；k++)
        Mh[k]＝0；
    wsum＝0；
    //遍历空间域邻域
    for(ii＝-hs；ii＜hs；ii++)
    {
        p＝i+ii；
        if(p＜0||p＞＝m_height)
            break；
        for(jj＝-hs；jj＜hs；jj++)
        {
            q＝j+jj；
            if(q＜0||q＞＝m_width)
                break；
            idxNei＝p*m_width+q；
            //特征向量
            ykNei[0]＝p；
            ykNei[1]＝q；
            for(k＝0；k＜m_dms；k++)
                ykNei[k+2]＝theInputData[idxNei+k*m_size]；
            //偏移向量
            for(k＝0；k＜m_dms+2；k++)
                shift[k]＝ykNei[k]-yk[k]；
        diff1＝fabs(shift[component1+2])；
        diff2＝fabs(shift[component2+2])；
        diff3＝fabs(shift[component3+2])；
        if(diff1＞hr||diff2＞hr||diff3＞hr)
            break；
        if(diff1＝＝0&&diff2＝＝0&&diff3＝＝0)
            break；
        w＝1；
        for(k＝0；k＜m_dms+2；k++)
            Mh[k]+＝w*shift[k]；
        wsum+＝w；
    }
}
if(wsum＞0)
{
    for(k＝0；k＜m_dms+2；k++)
    {
        Mh[k]＝Mh[k]/wsum；
        yk[k]+＝Mh[k]；
    }
}
else
{
    for(k＝0；k＜m_dms+2；k++)
    {
        Mh[k]＝0；
        flg＝false；
    }
                }
                disShift＝0；
                disShift+＝Mh[component1+2]；
                disShift+＝Mh[component2+2]；
                disShift+＝Mh[component3+2]；
                if(disShift＜3)
                    flg＝false；
                if(++nIter＝＝maxIter)
                    flg＝false；
            }
            for(k＝0；k＜m_dms；k++)
                theOutputData[idx+k*m_size]＝yk[k+2]；
         }
     }
    void SpatialclusterBasedSegmentAlgorithm::segment(double*theInputData，
double*theOutputData，long*cntData，int cntAttributes，long*labels，double
minDis)
    {
        long idx；
        int i，j，ii，jj，p，q；
        int k；
        long idxObj；
        long idxFeed，idxNei；
        long*cntPixels；
        int*indexTable；
        int index；
bool flg；
double delta，dis；
long cntObjects；
//初始化
cntPixels＝new long[m_size]；
indexTable＝new int[m_size]；
for(idx＝0；idx＜m_size；idx++)
{
    labels[idx]＝-1；
    for(k＝0；k＜cntAttributes；k++)
        theOutputData[idx+k*m_size]＝0；
}
//标号，填充对象集
idxObj＝-1；
for(idx＝0；idx＜m_size；idx++)
{
    //如果当前像元尚未标记
    if(labels[idx]＜0)
    {
        //用新的对象编号标记当前像元
        labels[idx]＝++idxObj；
        //当前对象包含的像元数为
        cntPixels[idxObj]＝1；
    }
    //如果该点已标记，则继续下一点
    else
        continue；
//以该像元作为种子点标记对象
idxFeed＝idx；
//索引表(已标记点编号-＞邻域点位置)
indexTable[0]＝idxFeed；
index＝0；
while(1)
{
    //行列号
    i＝idxFeed/m_width；
    j＝idxFeed％m_width；
    //没找到可标记的邻域点
    flg＝false；
    //遍历八邻域
    //ii和jj表示当前像元的相对位置，即与中心像元的偏移量
    //p和q表示当前像元的绝对位置，即在图像中的行列号
    for(ii＝-1；ii＜＝1；ii++)
    {
        for(jj＝-1；jj＜＝1；jj++)
        {
            //获取当前像元的绝对位置
            p＝i+ii；
            q＝j+jj；
            //如果超出图像边界，则继续
            if(p＜0||p＞＝m_height||q＜0||q＞＝m_width)
                 continue；
            //邻域像元编号
            idxNei＝p*m_width+q；
            //如果邻域点未被标记
            if(labels[idxNei]＝＝-1)
            {
               dis＝0；
                //如果邻域点与种子点的特征值相差较大则结束判
断
                for(k＝0；k＜cntAttributes；k++)
                {
                    delta＝theInputData[idxNei+k*m_size]-
theInputData[idxFeed+k*m_size]；
                    if(fabs(delta)＞minDis)break；
                }
                //特征向量相近
                if(k＝＝cntAttributes)
                {
                    //标记邻域点
                    labels[idxNei]＝idxObj；
                    //当前对象包含的像元数增加
                    cntPixels[idxObj]++；
                    //将邻域点放入索引表
                    indexTable[++index]＝idxNei；
                    //已找到可标记的邻域点
                    flg＝true；
                }
            }
        }
    }
//通过索引表进行深入搜索
//如果已找到可标记的邻域点，则以最后一个标记点为种子点继续搜索
//如果没找到可标记的邻域点，则以上一个标记点为种子点继续搜索
//如果已经没有上一个标记点，则区域填充完毕
        if(flg)
            idxFeed＝indexTable[index]；
        else if(index＞1)
            idxFeed＝indexTable[--index]；
        else
            break；//循环结束，完成填充
    }
}
//对象总数
cntObjects＝idxObj+1；
//赋值
for(idx＝0；idx＜m_size；idx++)
{
    //当前像元所属对象的编号
    idxObj＝labels[idx]；
            //遍历特征，累加
            for(k＝0；k＜cntAttributes；k++)
                 theOutputData[idxObj+k*cntObjects]+＝
theInputData[idx+k*m_size]；
        }
        //均值
        for(idxObj＝0；idxObj＜cntObjects；idxObj++)
        {
           for(k＝0；k＜cntAttributes；k++)
              theOutputData[idxObj+k*cntObjects]/＝cntPixels[idxObj]；
        }
        //对象数
        *cntData＝cntObjects；
    }
步骤4：对经步骤3平滑后的Landsat TM影像反射率图像中的影像第5波段进行OTSU直方图阈值分割，区分非高潮高地和高潮高地。实现方法的核心代码为：
    void CoralReefExtractAlgorithm::extractHighTideElevation(double*
theImageData，int*high)
    {
        //定义变量
        //分配内存
        mir＝new double[m_cntData]；
        mask＝new double[m_cntData]；
        //构建掩膜
        haveHigh＝false；
        cntMask＝0；
    for(idx＝0；idx＜m_cntData；idx++)
    {
        //获取中红外波段
        mir0＝theImageData[idx+4*m_cntData]；
        mir[idx]＝mir0；
        //判断是否存在高潮高地
        if(mir0＞150)
        {
             haveHigh＝true；
        }
        mask[cntMask]＝mir[idx]；
        cntMask++；
    }
    //计算阈值
    HistogramAssistProcessingTool*hap＝new
HistogramAssistProcessingTool()；
    hap-＞buildHistogram(mask，cntMask，100)；
    hap-＞soomth(2)；
    threshold＝hap-＞calThresholdByBreak()；
    //如果计算出的阈值过大，则掩膜后用大津法计算阈值
    if(threshold＜100)
    {
        cntMask＝0；
        for(idx＝0；idx＜m_cntData；idx++)
        {
            if(mir[idx]＞＝100)
        {
            mask[cntMask]＝mir[idx]；
            cntMask++；
        }
        }
        hap-＞buildHistogram(mask，cntMask，100)；
        threshold＝hap-＞calThresholdByOTSU()；
    }
    //分割提取
    for(idx＝0；idx＜m_cntData；idx++)
    {
        if(haveHigh&&mir[idx]＞threshold)
            high[idx]＝1；
        else
            high[idx]＝0；
    }
}
步骤5：对步骤4中的高潮高地，利用NDVI进行OTSU直方图阈值分割，区分陆地植被和灰沙/建筑。实现方法的核心代码为：
    void  CoralReefExtractAlgorithm::extractLandVegetationAndSandBuilding(double*
theImageData，int*high，int*landVegetation，int*sandOrBuilding)
    {
        //定义变量
        //分配内存
        ndvi＝new double[m_cntData]；
        mask＝new double[m_cntData]；
        //构建掩膜
       haveLandVegetation＝false；
        cntMask＝0；
        for(idx＝0；idx＜m_cntData；idx++)
        {
            //获取NDVI
        ndvi0＝(theImageData[idx+3*m_cntData]-theImageData[idx+2*m_cntData])
/(theImageData[idx+3*m_cntData]+theImageData[idx+2*m_cntData])；
        ndvi[idx]＝ndvi0；
        //判断是否存在陆地植被
        if(ndvi0＞0.4&&high[idx]＝＝1)
            haveLandVegetation＝true；
        if(high[idx]＝＝1)
        {
            mask[cntMask]＝ndvi[idx]；
            cntMask++；
        }
    }
    //计算阈值
    HistogramAssistProcessingTool*hap＝new
HistogramAssistProcessingTool()；
    hap-＞buildHistogram(mask，cntMask，100)；
    threshold＝hap-＞calThresholdByOTSU()；
    //分割提取
    for(idx＝0；idx＜m_cntData；idx++)
    {
        if(high[idx]＝＝1)
        {
            if(haveLandVegetation&&ndvi[idx]＞threshold)
            {
            landVegetation[idx]＝1；
            sandOrBuilding[idx]＝0；
        }
        else
        {
            landVegetation[idx]＝0；
            sandOrBuilding[idx]＝1；
        }
    }
    else
        landVegetation[idx]＝sandOrBuilding[idx]＝0；
    }
}
步骤6：对步骤4中的非高潮高地，利用影像第3波段进行直方图阈值分割，区分低潮高地和非干出质底。实现方法的核心代码为：
    void CoralReefExtractAlgorithm::extractLowTideElevation(double*
theImageData，int*high，int*low)
    {
        //定义变量
        //分配内存
        red＝new double[m_cntData]；
        mask＝new double[m_cntData]；
        //构建掩膜
        haveLow＝false；
        cntMask＝0；
        for(idx＝0；idx＜m_cntData；idx++)
       {
        //获取红光波段
        red0＝theImageData[idx+2*m_cntData]；
        red[idx]＝red0；
        //判断是否存在低潮高地
        if(red0＞100&&high[idx]＝＝0)
            haveLow＝true；
        if(high[idx]＝＝0)
        {
            mask[cntMask]＝red[idx]；
            cntMask++；
        }
    }
    //用断点法计算阈值
    HistogramAssistProcessingTool*hap＝new
HistogramAssistProcessingTool()；
    hap-＞buildHistogram(mask，cntMask，100)；
    hap-＞soomth(2)；
    threshold＝hap-＞calThresholdByBreak()；
    //如果计算出的阈值过大，则掩膜后用大津法计算阈值
    if(threshold＞150)
    {
        cntMask＝0；
        for(idx＝0；idx＜m_cntData；idx++)
        {
            if(high[idx]＝＝0&&red[idx]＜＝150)
            {
            mask[cntMask]＝red[idx]；
            cntMask++；
        }
    }
    hap-＞buildHistogram(mask，cntMask，100)；
    threshold＝hap-＞calThresholdByOTSU()；
}
//分割提取
for(idx＝0；idx＜m_cntData；idx++)
{
    if(haveLow&&high[idx]＝＝0&&red[idx]＞threshold)
        low[idx]＝1；
    else
        low[idx]＝0；
    }
}
步骤7：对步骤6中的低潮高地，利用NDVI进行直方图阈值分割，区分海藻/海草和浅水珊瑚。实现方法的核心代码为：
    void CoralReefExtractAlgorithm::extractAlgaeGrassAndShallowCorai(double*
theImageData，int*low，int*algaeOrGrass，int*shallowCoral)
    {
        //定义变量
        //分配内存
        ndvi＝new double[m_cntData]；
        nir＝new double[m_cntData]；
        mask＝new double[m_cntData]；
        //构建掩膜
        haveAlgaeOrGrass＝false；
    cntMask＝0；
    for(idx＝0；idx＜m_cntData；idx++)
    {
        //获取NDVI和近红外波段
    ndvi0＝(theImageData[idx+3*m_cntData]-theImageData[idx+2*m_cntData])
/(theImageData[idx+3*m_cntData]+theImageData[idx+2*m_cntData])；
        nir0＝theImageData[idx+3*m_cntData]；
        ndvi[idx]＝ndvi0；
        nir[idx]＝nir0；
        //判断是否存在海藻/海草
        if(ndvi0＞0.1&&nir0＞150&&low[idx]＝＝1)
              haveAlgaeOrGrass＝true；
        if(low[idx]＝＝1)
        {
            mask[cntMask]＝ndvi[idx]；
            cntMask++；
        }
    }
    //计算阈值
    HistogramAssistProcessingTool*hap＝new
HistogramAssistProcessingTool()；
    hap-＞buildHistogram(mask，cntMask，100)；
    threshold＝hap-＞calThresholdByOTSU()；
    //分割提取
    for(idx＝0；idx＜m_cntData；idx++)
    {
        if(low[idx]＝＝1)
       {
        if(haveAlgaeOrGrass&&ndvi[idx]＞threshold&&nir[idx]＞150)
        {
           algaeOrGrass[idx]＝1；
           shallowCoral[idx]＝0；
        }
        else
        {
           algaeOrGrass[idx]＝0；
           shallowCoral[idx]＝1；
        }
    }
    else
        algaeOrGrass[idx]＝shallowCoral[idx]＝0；
  }
}
步骤8：对步骤6中的非干出质底，利用影像第1波段进行直方图阈值分割，区分深水珊瑚和开放水体。实现方法的核心代码为：
    void CoralReefExtractAlgorithm::extractDeepCoralAndOpenWater(double*
theImageData，int*none，int*deepCoral，int*openWater)
    {
        //定义变量
        //分配内存
        blue＝new double[m_cntData]；
        mask＝new double[m_cntData]；
        //构建掩膜
        //haveAlgeaOrGrass＝false；
        cntMask＝0；
        for(idx＝0；idx＜m_cntData；idx++)
    {
        //获取NDVI
        blue0＝theImageData[idx+0*m_cntData]；
        blue[idx]＝blue0；
        if(none[idx]＝＝1)
        {
           mask[cntMask]＝blue[idx]；
           cntMask++；
        }
      }
    //计算阈值
    HistogramAssistProcessingTool* hap＝new
HistogramAssistProcessingTool()；
    hap-＞buildHistogram(mask，cntMask，100)；
    threshold＝hap-＞calThresholdByOTSU()；
    //分割提取
    for(idx＝0；idx＜m_cntData；idx++)
    {
       if(none[idx]＝＝1)
       {
           if(blue[idx]＞threshold)
           {
               deepCoral[idx]＝1；
               openWater[idx]＝0；
           }
           else
           {
               deepCoral[idx]＝0；
           openWater[idx]＝1；
         }
      }
      else
          deepCoral[idx]＝openWater[idx]＝0；
   }
}
在自主研发的海岸海洋遥感信息处理类库的支持下，采用标准C++实现了上述方法，应用该方法对八个试验区进行了珊瑚岛礁遥感信息自动提取，基本能够提取出陆地植被、灰沙/建筑、海藻/海草、浅水珊瑚、深水珊瑚5种珊瑚岛礁覆盖类型及开放水体背景信息，整体识别率优于90％。通过进一步分析得出，未能准确提取的信息主要包括：较厚的云团被误提取为灰沙/建筑，其原因是二者的光谱特征过于相近造成难以区分；较薄的云雾被误提取为海藻/海草、浅水珊瑚或深水珊瑚，其原因是较薄云雾的光谱中混合了云雾及其下垫面的信息，光谱成分复杂造成难以识别；云的阴影被误提取为开放水体，其原因同样为光谱特征相近；人工建筑被误提取为岛屿植被，其原因为岛屿上植被茂密，而人工建筑往往面积较小，因而混合像元情况比较严重。总体上，该方法经过少量后处理，能够应用于实际工作。

Claims (1)

1.一种基于Landsat数据源的珊瑚岛礁遥感信息自动提取方法，包括以下步骤：
步骤1：对原始遥感影像进行辐射定标，将DN值转为辐射亮度值，定标公式为L＝DN×gain+bias，增益系数gain和偏移量bias采用USGS提供的同期参数；
步骤2：对步骤1中定标后的遥感影像进行大气校正，将辐射亮度值转化为地表反射率，校正公式为：
R = ( L - L p ) · π ( E 0 · cos θ z · T z + E d ) · T v
其中，R为地表反射率，L为步骤1得出的辐射亮度值，E0为大气顶端的太阳辐照度，可由通过太阳常数和日地距离求得，即E0=Es/D2，θz为太阳天顶角，LP为大气的程辐射能量值，Tz为从太阳到地表入射方向上的大气透射率，Tv为从地表到传感器反射方向上的大气透射率，Ed为经散射作用到达地表的太阳辐照度，具体参数Tz、Tv、Ed、LP的计算公式为：
Tz＝cosθz
Tv＝cosθv
Ed＝πLp
Lp＝L-0.02·(E0·cosθz·Tz+Ed)·Tv
其中，θv为传感器的天顶角；
步骤3：对经步骤2大气校正后的影像进行Mean Shift空间滤波，以减轻条带和椒盐噪声，得到平滑后的TM反射率图像；
步骤4：对经步骤3平滑后的TM反射率图像中的TM5波段进行Otsu直方图阈值分割，区分非高潮高地和高潮高地；
步骤5：对步骤4中的高潮高地，利用NDVI进行Otsu直方图阈值分割，区分陆地植被和灰沙/建筑；
步骤6：对步骤4中的非高潮高地，利用TM3波段进行直方图阈值分割，区分低潮高地和非干出质底；
步骤7：对步骤6中的低潮高地，利用NDVI进行直方图阈值分割，区分海藻/海草和浅水珊瑚；
步骤8：对步骤6中的非干出质底，利用TM1波段进行直方图阈值分割，区分深水珊瑚和开放水体。
CN2010105675077A 2010-12-01 2010-12-01 基于Landsat数据源的珊瑚岛礁遥感信息自动提取方法 Expired - Fee Related CN102032903B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN2010105675077A CN102032903B (zh) 2010-12-01 2010-12-01 基于Landsat数据源的珊瑚岛礁遥感信息自动提取方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN2010105675077A CN102032903B (zh) 2010-12-01 2010-12-01 基于Landsat数据源的珊瑚岛礁遥感信息自动提取方法

Publications (2)

Publication Number Publication Date
CN102032903A CN102032903A (zh) 2011-04-27
CN102032903B true CN102032903B (zh) 2012-11-21

Family

ID=43886105

Family Applications (1)

Application Number Title Priority Date Filing Date
CN2010105675077A Expired - Fee Related CN102032903B (zh) 2010-12-01 2010-12-01 基于Landsat数据源的珊瑚岛礁遥感信息自动提取方法

Country Status (1)

Country Link
CN (1) CN102032903B (zh)

Families Citing this family (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102231190B (zh) * 2011-07-08 2012-10-31 中铁第四勘察设计院集团有限公司 冲洪积扇信息的自动提取方法
CN102750683B (zh) * 2012-06-18 2014-10-29 常州大学 Modis遥感影像中海面条带噪声和条状云的过滤方法
CN103177255B (zh) * 2013-03-15 2016-01-20 浙江大学 一种基于多分辨率数字高程模型的潮间带提取方法
CN104318544B (zh) * 2014-09-25 2017-08-11 中国水产科学研究院东海水产研究所 基于夜间灯光卫星遥感数据光诱捕作业渔船数量估算方法
CN104794724A (zh) * 2015-05-04 2015-07-22 福建师范大学 一种基于非线性尺度滤波的遥感影像建筑物提取方法
CN104820972B (zh) * 2015-05-07 2017-06-16 北京空间机电研究所 一种基于在轨分类统计的红外影像me噪声去除方法
CN109635765B (zh) * 2018-12-19 2023-08-29 海南空天信息研究院 一种浅海珊瑚礁遥感信息自动提取方法
CN110411962B (zh) * 2019-08-03 2022-04-08 国家海洋环境监测中心 珊瑚礁白化遥感监测机理分析方法
CN110544236B (zh) * 2019-08-03 2023-06-16 国家海洋环境监测中心 基于时间序列卫星影像的珊瑚礁白化遥感监测方法
CN111695503B (zh) * 2020-06-11 2023-04-18 自然资源部第一海洋研究所 基于双波段辐亮度的珊瑚礁底质分类方法
CN115078263B (zh) * 2022-05-27 2023-07-28 苏州科技大学 顾及潮汐影响的海草遥感信息提取方法

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1854757A (zh) * 2005-04-28 2006-11-01 中国科学院遥感应用研究所 遥感影像群判读系统方法
CN101034472A (zh) * 2007-03-29 2007-09-12 上海大学 Gis支持下的卫星遥感数字图像的地形变换
CN101871884A (zh) * 2010-06-02 2010-10-27 中国国土资源航空物探遥感中心 多景aster遥感数据大气校正与区域性矿物填图方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN1854757A (zh) * 2005-04-28 2006-11-01 中国科学院遥感应用研究所 遥感影像群判读系统方法
CN101034472A (zh) * 2007-03-29 2007-09-12 上海大学 Gis支持下的卫星遥感数字图像的地形变换
CN101871884A (zh) * 2010-06-02 2010-10-27 中国国土资源航空物探遥感中心 多景aster遥感数据大气校正与区域性矿物填图方法

Non-Patent Citations (8)

* Cited by examiner, † Cited by third party
Title
SPOT、QuickBird卫星遥感数据提取东沙岛植被信息的比较;初佳兰等;《海洋学研究》;20060630(第02期);79-85 *
冷秀华等.海岛遥感图像数据管理示范系统的构架与分析.《海岸工程》.2007,(第03期),72-77.
初佳兰等.SPOT、QuickBird卫星遥感数据提取东沙岛植被信息的比较.《海洋学研究》.2006,(第02期),79-85.
基于光学遥感的海岛潮间带和湿地信息提取――以东沙岛(礁)为例;王小龙等;《海洋科学进展》;20051015(第04期);477-481 *
海岛遥感图像数据管理示范系统的构架与分析;冷秀华等;《海岸工程》;20070915(第03期);72-77 *
王小龙等.基于光学遥感的海岛潮间带和湿地信息提取――以东沙岛(礁)为例.《海洋科学进展》.2005,(第04期),477-481.
遥感技术应用于海岛保护与利用规划的可行性研究;马毅等;《海洋开发与管理》;20090715(第07期);92-95 *
马毅等.遥感技术应用于海岛保护与利用规划的可行性研究.《海洋开发与管理》.2009,(第07期),92-95.

Also Published As

Publication number Publication date
CN102032903A (zh) 2011-04-27

Similar Documents

Publication Publication Date Title
CN102032903B (zh) 基于Landsat数据源的珊瑚岛礁遥感信息自动提取方法
Schulz et al. Land use mapping using Sentinel-1 and Sentinel-2 time series in a heterogeneous landscape in Niger, Sahel
CN109190538A (zh) 一种基于遥感技术的多泥沙河流三角洲海岸带演化分析方法
Zhang et al. GWL_FCS30: global 30 m wetland map with fine classification system using multi-sourced and time-series remote sensing imagery in 2020
Baumstark et al. Mapping seagrass and colonized hard bottom in Springs Coast, Florida using WorldView-2 satellite imagery
CN103063311B (zh) 基于土壤指数的裸露基岩信息提取方法
CN103646246A (zh) 一种基于决策树模型的多光谱遥感影像河流信息提取方法
CN102855494A (zh) 一种卫星遥感影像的水体提取方法及装置
Lu et al. Detection of urban expansion in an urban-rural landscape with multitemporal QuickBird images
CN113971769B (zh) 基于多源大数据的海岸带地域功能长时序识别方法
St-Pierre et al. Kelp-bed dynamics across scales: Enhancing mapping capability with remote sensing and GIS
Li Textural and rule-based lithological classification of remote sensing data, and geological mapping in Southwestern Prieska sub-basin, Transvaal Supergroup, South Africa
CN113780307A (zh) 一种区域年度最大蓝绿空间信息提取方法
Whiteside et al. A semi-automated approach for quantitative mapping of woody cover from historical time series aerial photography and satellite imagery
Suardana et al. Estimation and Mapping Above-Ground Mangrove Carbon Stock Using Sentinel-2 Data Derived Vegetation Indices in Benoa Bay of Bali Province, Indonesia
Wang et al. Extraction of palaeochannel information from remote sensing imagery in the east of Chaohu Lake, China
Li Dynamic monitoring algorithm of natural resources in scenic spots based on MODIS Remote Sensing technology
Civco et al. A hierarchical approach to land use and land cover mapping using multiple image types
Lan et al. Spectral index-spatially correlated water extraction method based on GF-5 hyperspectral satellite images
Afrasinei Study of land degradation and desertification dynamics in North Africa areas using remote sensing techniques
CN115019184B (zh) 基于遥感影像的石漠化程度自动分级方法及装置
Aravena Satellite Imagery Interpretation for Environmental Mapping and Monitoring
Hong et al. Measuring canopy morphology of saltmarsh plant patches using UAV-based LiDAR data
CN115049947A (zh) 一种基于无人机多光谱数据的河床基岩识别方法、系统、设备及应用
CN116797921A (zh) 基于Sentinel-2影像的山区湖泊水体提取方法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
C14 Grant of patent or utility model
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20121121

Termination date: 20151201

EXPY Termination of patent right or utility model