CN105259145A - 一种同时遥感岛礁水下地形和地物的方法 - Google Patents

Abstract

本发明公开了一种同时遥感岛礁水下地形和地物的方法，包括以下步骤：(A)建立水面遥感反射率与水深、水下地物反射率的物理分析模型；(B)测量水体吸收系数、水体散射系数以及研究区典型地物的反射率；(C)对卫星原始图像进行大气纠正，并消除水面镜面反射光的影响；(D)根据所述物理分析模型和测量得到的水体吸收系数、水体散射系数和某一典型地物反射率，计算得到某一像元处的水深、水下地物类型；(E)重复步骤(D)，得到所有像元的水深、水下地物类型。该方法可以实现水深在40米以内的清洁浅水水域水深和水下地物类型的同时遥感反演，且具有较高的精度。可对浅海和岛礁水下地形、地物及生态环境进行遥监测。

Description

一种同时遥感岛礁水下地形和地物的方法

技术领域

本发明属于卫星遥感技术领域，特别涉及一种同时遥感岛礁水下地形和地物的方法。

背景技术

采用卫星遥感数据进行岛礁浅水区水下地形和地物调查具有效率高，不受水面条件限制等优点。如方法得当，还可以较客观细致地反映浅水区水下地形和地物分布，对海岛和浅海区的经济建设、生态调查，乃至军事应用均具有良好前景。

目前采用遥感数据进行水下地形地物调查有三类方法：

第一类，是采用经验公式方法。即实地测量不同水深和不同地物的水面光谱，建立水深、水下地物类型与水面反射率的经验公式或数据库，再根据遥感数据，用该经验公式估算水下地形地物，或对数据库用查表或光谱匹配方法估算水下地形地物。

但是，采用经验公式的方法建立的模型物理机理不清，需要进行大量的水面与水下地物调查，且适应性不好。最为严重的是，成像条件千差万别，事先很难完全测量。故这类方法精度不高。

第二类，是半经验方法。这类方法总体上根据水体辐射传输原理建模，但对于模型的难以获取的部分，如成像时的环境条件因子，水的消光系数等，用经验的方法求得。

采用半经验、半分析模型的方法虽然建立的模型具有一定的物理意义，但也在一定程度上具有经验公式的不足，方法只适用于特定的情况。由于很难用经验方法模拟所有的实际情况，这类方法精度和适应性也受到相当的限制。

第三类，是建立物理分析模型的方法。这种方法根据水体辐射传输机理和过程建模，有清晰的物理学意义，如测量了水和杂质的吸收、散射系数和水下地物反射率等稳定的物理参数，将具有较好的精度和适应性。是该领域研究的发展方向。

然而，对于目前的水下地形和地物遥感，即使建立了正确的水体反射率与水深和水下地物的物理分析模型，也存在以下方面问题。

其一，对于水下地形和地物遥感，影响最大的成像条件因素为大气和水面镜面反射率影响。大气的影响通过大气纠正去除，这方面已进行了大量的研究，已有较为通用有效的方法；水面镜面反射的影响目前只通过选定太阳入射和观测方向避开严重的太阳耀斑，对一般的水面镜面反射则无法消除，故有人采用经验公式的方法处理该问题。然而，镜面反射的强度与水面粗糙度有关，每个点皆有不同，故这种方法无法取得理想的效果。

其二，水面反射率取决于水深和水下地物反射率，由于对一般的地物而言，每个波段的反射率皆不相同，加上水深，故不管采用多少个波段的遥感数据，皆是未知数多于方程的个数，无法求解。故现有的方法多假定水下反射率已知，也即只适用于特定的水域。

其三，不同的遥感数据波段的透水能力不同，多是采用透水能力较强的兰或绿波段来进行水深反演。由于透水能力强，在浅水区的水深不敏感，水深的精度不高；有的采用红波段或红外波段(或其比值)来进行反演，只能反演较浅的水深。对水深较大，如水深大于5米时，则无法进行水深和水下地物反演。

发明内容

本发明的目的是提供一种同时遥感岛礁水下地形和地物的方法，该方法可以实现水深在40米以内的清洁浅水水域水深和水下地物类型的同时遥感反演，且具有较高的精度。可对浅海和岛礁水下地形、地物及生态环境进行遥监测。

为了达到上述目的，本发明采用如下技术方案：

一种同时遥感岛礁水下地形和地物的方法，包括以下步骤：

(A)建立水面遥感反射率与水深、水下地物反射率的物理分析模型。

(B)测量水体吸收系数、水体散射系数以及研究区典型地物的反射率。

(C)对卫星原始图像进行大气纠正，并消除水面镜面反射光的影响。

(D)对某一像元根据所述物理分析模型和步骤(B)中测量得到的水体吸收系数、水体散射系数和某一典型地物反射率计算得到该像元处的水深、水下地物类型。

(E)重复步骤(D)，得到得到所有像元的水深、水下地物类型。

作为一种具体的实施例，所述步骤(A)中水面遥感反射率与水深、水下地物反射率的物理分析模型的函数关系为：

$R_w=\frac{\mathrm{\β}P\left(\mathrm{\Θ}\right)}{4\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})}(1-e^{-\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})h})+R_b{e}^{-\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})h};$

式中：P(Θ)为水体散射相函数，Θ为水体散射角_；

β为水体散射系数，α为水体吸收系数，R_b为水底底质反射率，h为水深；

μ＝secθ+secφ，θ为阳光在水中传播的天顶角，φ为水中的观测天顶角；

R_w＝R_ws+R_wb，R_w为水体反射率，R_ws为水体散射光对反射率的贡献，R_wb为水底反射光对反射率的贡献。

作为一种具体的实施例，根据权利要求1所述的同时遥感岛礁水下地形和地物的方法，其特征在于，步骤(C)具体为：

(C1)从图像上提取暗像元，迭代计算大气散射系数、大气透过率和天空光辐照度。

(C2)对图像中所有像元进行大气纠正。

(C3)消除天空光镜面反射光。

(C4)消除太阳直射光镜面反射光。

进一步地，所述步骤(C1)中迭代计算大气散射系数、大气透过率和天空光辐照度的步骤包括：

(C11)读取暗像元的表观反射率R_vd，根据暗像元的表观反射率R_vd的计算公式：

$R_{vd}=\frac{{\mathrm{\ωTR}}_v}{2}+\frac{\mathrm{\ω}}{4}P\left(\mathrm{\θ}\right)---\left(1\right)$

式中：ω为大气散射系数，T为大气透过率，R_v为地物反射率，P(θ)为大气散射相函数，θ为散射角，

首先忽略等式(1)右边第一项的值，计算得到大气散射系数的初值

(C12)：根据大气透过率的初值的计算公式：T′＝e^-ω′，计算出大气透过率的初值T′，将大气透过率的初值T′代入等式(1)右边第一项中的T，将大气散射系数的初值ω′代入等式(1)右边第一项中的ω，计算得到大气散射系数

(C13)：计算天空光幅照度：

将步骤(C12)中计算得到的大气散射系数ω₁代入天空光幅照度计算公式计算天空光幅照度E_sky1；

式中：c为大气纠正前后的大气散射系数比，为测量值；

(C14)：将步骤(C12)中计算得到的大气散射系数ω₁代入大气透过率的计算公式T＝e^-ω，计算出大气透过率T₁；

(C15)：将步骤(C12)中计算得到的大气散射系数ω₁代入等式(1)右边第一项中的ω、将步骤(C14)中计算得到的大气透过率T₁代入等式(1)右边第一项中的T，计算出精度高于大气散射系数ω₁的大气散射系数ω₂；

根据步骤(C13)计算出精度更高的天空光幅照度E_sky2；

根据步骤(C14)计算出精度高于大气透过率T₁的大气透过率T₂；

依次办法，迭代计算出精度高于大气散射系数ω₂的大气散射系数ω_k、精度高于大气透过率T₂的大气透过率T_k和精度高于天空光幅照度E_sky2的天空光幅照度E_skyk；其中，k>2。

进一步地，所述步骤(C2)中对图像中所有像元进行大气纠正的计算公式为：

$R_g=\frac{R_1-\frac{{\mathrm{\ω}}_k}{4}P\left(\mathrm{\θ}\right)}{\[T^{\′\′}+\frac{{\mathrm{\ω}}_k}{2}\]T_k};$

式中：R_g为大气纠正后得到的像元反射率，T″＝e^-τ/cosθ,其中τ为大气光学厚度，R₁为像元表观反射率：L为像元原始辐亮度，E₀为太阳辐照度。

进一步地，所述步骤(C3)的具体步骤为：将像元各波段的原始辐亮度减去步骤(A5)得到的天空光辐照度E_skyk与天空光水面镜面反射率R_m之积，其公式为：L₁＝L-E_skykR_m，

式中：L₁为像元消除天空光镜面反射后的辐亮度。

进一步地，所述步骤(C4)的具体步骤为：将图像转化为反射率图像，像元所有波段的反射率均减去红外波段的反射率，其中，像元所有的波段的反射率减去红外波段的反射率的公式为：R₂＝R_g-R_ir

式中：R₂为像元消除光镜面反射后的反射率，R_ir为红外波段的反射率。

作为一种具体的实施例，所述步骤(D)中水深的计算公式为：

$h=\frac{1}{\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})}ln\[\frac{4\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})R_b-\mathrm{\β}P\left(\mathrm{\Θ}\right)}{4\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})R_w-\mathrm{\β}P\left(\mathrm{\Θ}\right)}\];$

式中：P(Θ)为水体散射相函数，Θ为水体散射角_；

β为水体散射系数，α为水体吸收系数，R_b为水底底质反射率，h为水深；

μ＝secθ+secφ，θ为阳光在水中传播的天顶角，φ为水中的观测天顶角；

R_w＝R_ws+R_wb，R_w为水体反射率，R_ws为水体散射光对反射率的贡献，R_wb为水底反射光对反射率的贡献。

作为一种具体的实施例，步骤(D)具体为：

(D1)将步骤(B)中测量得到的某一典型地物反射率、所有波段的遥感数据代入所述物理分析模型，计算得到若干不同的水深，计算得到水深初值；根据计算的水深初值和水深计算最佳波段规则，确定水深计算最佳波段，根据水深计算最佳波段的遥感数据算出最优水深值。

所述水深计算最佳波段规则为：当水深值在(0，0.3)米时，最佳波段为近红外波段；当水深值在[0.3，5)米时，最佳波段为近红光波段；当水深值在[5，20)米时，最佳波段为绿光波段；当水深值在[20，45]米时，最佳波段为兰光波段。

(D2)将所有波段的遥感数据、所述最优水深值代入所述物理分析模型，计算得到各波段的地物反射率，然后将计算得到的各波段的地物反射率与测量得到的地物反射率进行匹配，计算其误差。

其中，具有最小误差的地物反射率所对应的地物为该像元对应的水下地物类型，根据该地物反射率计算得到最佳水深为该像元的水深。

进一步地，所述步骤(D2)中计算得到的各波段的地物反射率与测量得到的地物反射率的误差的计算公式为：

${\mathrm{\ϵ}}^2=\underset{i=n}{\Σ}{(R_b-R_b^{\′})}^2;$

式中：R_b为测量得到的地物反射率，R′_b为计算得到的波段地物反射率。

与现有技术相比，本发明的有益效果在于：

本发明一种同时遥感岛礁水下地形和地物的方法，首先通过，实现水深在40米以内的清洁浅水水域水深和水下地物类型的同时遥感反演，且结果具有较高的精度。可对浅海和岛礁水下地形、地物及生态环境进行遥监测。

通过对西沙永兴岛水下地形和地物遥感反演，很好地验证了本发明方法在进行水下地形和地物反演的过程中具有较高的精度。

附图说明

图1是本发明的同时遥感岛礁水下地形和地物的方法流程图。

图2是本发明的岛礁浅水区水深与水下地物遥感模型的辐射传输原理示意图。

图3是400nm到2500nm海水吸收光谱曲线图。

图4是海水反射率曲线图。

图5a是原始卫星图像。

图5b是经大气纠正后的卫星图像。

图5c是消除水面镜面反射光后的卫星图像。

图6是永兴岛珊瑚礁海域水深遥感图。

图7a、7b为永兴岛珊瑚礁海域水深遥感图三维立体图。

图8是永兴岛水下地物分布图。

图9是永兴岛珊瑚礁海域水深遥感图与普通海图的对比图。

图10是采用双通道测深仪对永兴岛进行水深测量，4条测线的分布图。

图11a—11d是图10中的4条测线的剖面的测量值与遥感值的对比曲线图。

具体实施方式

下面将结合附图以及具体实施方法来详细说明本发明，在本发明的示意性实施及说明用来解释本发明，但并不作为对本发明的限定。

如图1所示，本发明一种同时遥感浅水区地形和地物的方法，包括以下步骤：

(A)：建立水面遥感反射率与水深、水下地物反射率的物理分析模型。

如图2所示，为本发明的岛礁浅水区水深与水下地物遥感模型的辐射传输原理示意图。卫星传感器接收到的水面光能量由三部分组成，其一为水面镜面反射光，其二为水体散射光，其三为水底反射光。

其中，水面镜面反射光没有水深和水下地物的信息，为干扰信息，通过消除水域遥感数据的水面镜面反射影响的方法，去镜面反射光予以消除。消除水面镜面反射后，水体反射率R_w可由下式表示：

R_w＝R_ws+R_wb——(1)；

式中：R_ws为水体散射光对反射率的贡献，R_wb为水底反射光对反射率的贡献。

其中，水体散射光对反射率的贡献R_ws的计算公式为：

$R_{ws}=\frac{\mathrm{\β}P\left(\mathrm{\Θ}\right)}{4\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})}(1-e^{-\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})h})---\left(2\right);$

水底反射光对反射率的贡献R_wb为：

R_wb＝R_be^-μ(α+β)h——(3)；

将(2)式、(3)式代入(1)式中，得水面遥感反射率与水深、水下地物反射率的函数关系为：

$R_w=\frac{\mathrm{\β}P\left(\mathrm{\Θ}\right)}{4\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})}(1-e^{-\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})h})+R_b{e}^{-\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})h}---\left(4\right);$

式中：P(Θ)为水体散射相函数，Θ为水体散射角_，

β为水体散射系数，α为水体吸收系数，R_b为水底底质反射率，h为水深；

μ＝secθ+secφ；θ为阳光在水中传播的天顶角，φ为水中的观测天顶角。

(4)式即为水面遥感反射率与水深、水下地物反射率的物理分析模型。

(B)：精确测量水体消光系数、吸收系数和反向散射系数，以及研究区所有典型地物的反射率。根据测量得到的水体消光系数、吸收系数和反向散射系数算出水体吸收系数、散射系数。所述典型地物包括有沙、礁石、石珊瑚、软珊瑚、死珊瑚、白化珊瑚、海草等。

(C)：对卫星原始图像进行精确的大气纠正并消除水面镜面反射光的影响。

(D)：将典型地物反射率代入所述水面遥感反射率与水深、水下地物反射率的物理分析模型，计算水深。

假设某一像元处的地物反射率已知，将测得的典型地物反射率代入所述物理分析模型，根据不同波段遥感数据计算，计算得到若干个不同的水深初值，根据计算的水深初值和水深计算最佳波段规则，确定水深计算最佳波段，根据水深计算最佳波段的遥感数据算出最优水深值。

所述水深计算最佳波段规则为：当水深值在(0，0.3)米时，最佳波段为近红外波段；当水深值在[0.3，5)米时，最佳波段为近红光波段；当水深值在[5，20)米时，最佳波段为绿光波段；当水深值在[20，45]米时，最佳波段为兰光波段。具体如下表所示：

水深(米)	<0.3	0.3-5	5-20	20-45
					最优波段	近红外	红光	绿光	兰光

将所述最优水深值代入所述物理分析模型，同时代入所有波段的遥感数据，计算得到各波段的地物反射率；

将得到的各波段的地物反射率与测量得到的地物反射率进行匹配，误差最小的即是地物反射率所对应的地物为该像元对应的水下地物类型，根据该地物反射率计算得到的最佳水深为该像元的水深。

具体的，对不同的水深范围，根据各波段的消光强度，用不同的波段来计算。

将典型地物反射率带入模型即相当于地物反射率已知，水深计算公式为：

$h=\frac{1}{\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})}ln\[\frac{4\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})R_b-\mathrm{\β}P\left(\mathrm{\Θ}\right)}{4\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})R_w-\mathrm{\β}P\left(\mathrm{\Θ}\right)}\];$

(E)：将所有波段的遥感数据带入，根据上步计算得到的最优水深，计算各波段地物的反射率：

$R_b=R_w{e}^{\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})h}-\frac{\mathrm{\β}P\left(\mathrm{\Θ}\right)}{4\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})}(e^{\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})h}-1);$

然后将计算的各波段反射率与测量得到的地物反射率进行匹配，计算其误差。

${\mathrm{\ϵ}}^2=\underset{i=n}{\Σ}{(R_b-R_b^{\′})}^2;$

R_b为地物反射率，R′_b为计算的波段反射率。

(F)：将具有最小ε²的水下地物类型和水深作为像元的水下地物类型和水深。对所有水域像元均作该计算，得到水下地形图和水下地物分布图。

其中，步骤(C)对卫星原始图像进行精确的大气纠正并消除水面镜面反射光的影响，包括以下步骤：

步骤(C1)：从图像上提取暗像元，迭代计算大气散射系数、大气透过率和天空光辐照度。

首先，从原始卫星图像上提取暗像元。根据提取的像元，计算获取大气散射系数、大气透过率和天空光辐照度。具体步骤如下：

(C11)读取暗像元的表观反射率R_vd，计算大气散射系数初值ω′：根据暗像元的表观反射率R_vd的计算公式： $R_{vd}=\frac{{\mathrm{\ωTR}}_v}{2}+\frac{\mathrm{\ω}}{4}P\left(\mathrm{\θ}\right)---\left(1\right)$

首先忽略等式(1)右边第一项的值，计算得到大气散射系数的初值

式中：ω为大气散射系数，T为大气透过率，R_v为地物反射率，P(θ)为大气散射相函数，θ为散射角。

(C12)：根据大气透过率的初值的计算公式：T′＝e^-ω′，计算出大气透过率的初值T′，并将大气透过率的初值T′代入等式(1)右边第一项中的T，将大气散射系数的初值ω′代入等式(1)右边第一项中的ω，计算得到大气散射系数

(C13)：计算天空光幅照度：

将步骤(C12)中计算得到的大气散射系数ω₁代入天空光幅照度计算公式计算天空光幅照度E_sky1，式中：c为大气纠正前后的大气散射系数比，通过测量得到，为小于2的测量值。

(C14)：将步骤(C12)中计算得到的大气散射系数ω₁代入大气透过率的计算公式T＝e^-ω，计算出大气透过率T₁。

(C15)：将步骤(C12)中计算得到的大气散射系数ω₁代入等式(1)右边第一项中的ω、将步骤(C14)中计算得到的大气透过率T₁代入等式(1)右边第一项中的T，计算出精度高于大气散射系数ω₁的大气散射系数ω₂。

根据步骤(C13)计算出精度更高的天空光幅照度E_sky2。

根据步骤(C14)计算出精度高于大气透过率T₁的大气透过率T₂。

依次办法，迭代计算出精度高于大气散射系数ω₂的大气散射系数ω_k、精度高于大气透过率T₂的大气透过率T_k和精度高于天空光幅照度E_sky2的天空光幅照度E_skyk；其中，k>2。

步骤(C2)：对图像中所有像元进行大气纠正。对图像中所有像元，均进行以下计算： $R_g=\frac{R_1-\frac{{\mathrm{\ω}}_k}{4}P\left(\mathrm{\θ}\right)}{\[T^{\′\′}+\frac{{\mathrm{\ω}}_k}{2}\]T_k},$

式中：R_g为大气纠正后得到的像元反射率，T′＝e^-τ/cosθ，其中τ为大气光学厚度，R₁为像元表观反射率：L为像元原始辐亮度，E₀为太阳辐照度。

经过大气纠正后，大气散射光的影响被很好消除，影像清晰度得到明显提高。但水域镜面反射光的影响严重，水质及水下地物信息受到严重干扰。

本发明根据水面镜面反射的特征，消除每个像元的水面镜反射，其中，水面镜反射由水面对太阳直射光的镜面反射和对天空光的镜面反射光组成。由于太阳直射光和天空光的光谱组成不同，故对二者的镜面反射光分别消除。

步骤(C3)消除天空光镜面反射光。

通过测量或暗像元提取等方法获取天空光强度，然后像元各波段辐亮度减去天空光辐照度与镜面反射率之积即可消除天空光镜面反射的影响。

将像元各波段的原始辐亮度减去步骤(A5)得到的天空光辐照度E_skyk与天空光水面镜面反射率R_m之积，其公式为：L₁＝L-E_skykR_m。

式中：L₁为像元消除天空光镜面反射后的辐亮度。

(C4)消除太阳直射光镜面反射光。

根据大量的测量实验证明，各光学波段的水面镜面反射光的反射率近相等。对同一像元，各波段的镜面反射率近相等；长波近红外波段，由于水强吸收，从水中出来的光强为0，其像元反射率即为镜面反射光的反射率。像元所有波段的反射率均减去红外波段的反射率，即可消除太阳直射光镜面反射的影响。

将图像转化为反射率图像，像元所有波段的反射率均减去红外波段的反射率，其中，波段的反射率减去红外波段的反射率的公式为：R₂＝R-R_ir。

式中：R像元原始反射率，R₂像元消除光镜面反射后的反射率，R_ir红外波段原始反射率。

根据大气纠正后得到的像元发射率R_g，则得到R₂＝R_g-R_ir。以此消除太阳直射光镜面反射的影响。

消除镜面反射后，图像水域部分的清晰度得到显著改善，镜面反射的影响得到有效消除，水下及水质信息得到清楚显示。效果良好。

采用原始数据中阴影和清深水体等暗像元计算天空光强度的方法，无需实地测量，可通过计算机自动完成水面镜面反射光的消除，更有实用价值。用叠代方法计算可得到较精确的天空光辐照度，从而更精确地消除天空光镜面反射的影响。

应用实例：

通过对该方法进行大量的应用试验，包括对西沙永兴岛水下地形和地物遥感反演，结果表明具有较高的精度。

1、测量西沙海域水体吸收系数、散射系数和典型地物反射率。部分测量结果如图3、图4所示。

2、进行精确的大气纠正并消除水面镜面反射光的影响

采用美国的WorldView2卫星遥感数据。如图5a所示，原始卫星图像大气散射严重，散射兰光明显影响了图像的清析度。如图5b所示，在经过大气纠正后，大气散射的影响已被告消除，清析度得到提高。但水域的水面镜面反射仍严重，水下地物信息受到严重干扰。如图5c所示，在经过大气纠正后，进一步消除镜面反射的影响后得到的图像，可看到水域的水面镜面反射已基本都被消除，水下于物信息清楚。

3、水下地形和水下地物的遥感反演

应用以上方法，计算得到每个像元的水深与水下地物类型，得到水下地形图和水下地物分布图如图6、图7a、图7b和图8所示。

4、精度验证

1)水下地形的反演精度验证。

(1)与海图的对比

如图9所示，对比遥感的水下地形图与海图。可见反演的地形图与海图具有较好的一致性，其地形的细节远比海图丰富。

(2)与实测水深的对比

采用双通道测深仪对永兴岛进行水深测量，共测量4条测线L1、L2、L3和L4，分布如图10所示。如图11a—11d所示，对比4条剖面的测量值与遥感值，水深测量值与遥感反演值非常接近，平均相对误差为18％。

2)水下地物的反演精度验证

选取35个样点，潜入水下调查水下地物类型。将遥感反演的地物类型与实际地物类型进行比较，仅有6个出现错判，准确率为83％。

以上对本发明实施例所提供的技术方案进行了详细介绍，本文中应用了具体个例对本发明实施例的原理以及实施方式进行了阐述，以上实施例的说明只适用于帮助理解本发明实施例的原理；同时，对于本领域的一般技术人员，依据本发明实施例，在具体实施方式以及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。

Claims (10)

1.一种同时遥感岛礁水下地形和地物的方法，其特征在于，包括以下步骤：

(A)建立水面遥感反射率与水深、水下地物反射率的物理分析模型；

(B)测量水体吸收系数、水体散射系数以及研究区典型地物的反射率；

(C)对卫星原始图像进行大气纠正，并消除水面镜面反射光的影响；

(D)根据所述物理分析模型和步骤(B)中测量得到的水体吸收系数、水体散射系数和某一典型地物反射率，计算得到某一像元处的水深、水下地物类型；

(E)重复步骤(D)，得到所有像元的水深、水下地物类型。

2.根据权利要求1所述的同时遥感岛礁水下地形和地物的方法，其特征在于，所述步骤(A)中水面遥感反射率与水深、水下地物反射率的物理分析模型的函数关系为：

$R_w=\frac{\mathrm{\β}P\left(\mathrm{\Θ}\right)}{4\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})}(1-e^{-\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})h})+R_b{e}^{-\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})h};$

式中：P(Θ)为水体散射相函数，Θ为水体散射角；

β为水体散射系数，α为水体吸收系数，R_b为水底底质反射率，h为水深；

μ＝secθ+secφ，θ为阳光在水中传播的天顶角，φ为水中的观测天顶角；

R_w＝R_ws+R_wb，R_w为水体反射率，R_ws为水体散射光对反射率的贡献，R_wb为水底反射光对反射率的贡献。

3.根据权利要求1所述的同时遥感岛礁水下地形和地物的方法，其特征在于，步骤(C)具体为：

(C1)从图像上提取暗像元，迭代计算大气散射系数、大气透过率和天空光辐照度；

(C2)对图像中所有像元进行大气纠正；

(C3)消除天空光镜面反射光；

(C4)消除太阳直射光镜面反射光。

4.根据权利要求3所述的同时遥感岛礁水下地形和地物的方法，其特征在于，所述步骤(C1)中迭代计算大气散射系数、大气透过率和天空光辐照度的步骤包括：

(C11)读取暗像元的表观反射率R_vd，根据暗像元的表观反射率R_vd的计算公式：

$R_{vd}=\frac{{\mathrm{\ωTR}}_v}{2}+\frac{\mathrm{\ω}}{4}P\left(\mathrm{\θ}\right)---\left(1\right)$

式中：ω为大气散射系数，T为大气透过率，R_v为地物反射率，P(θ)为大气散射相函数，θ为散射角，

首先忽略等式(1)右边第一项的值，计算得到大气散射系数的初值

(C12)：根据大气透过率的初值的计算公式：T′＝e^-ω′，计算出大气透过率的初值T′，将大气透过率的初值T′代入等式(1)右边第一项中的T，将大气散射系数的初值ω′代入等式(1)右边第一项中的ω，计算得到大气散射系数

(C13)：计算天空光幅照度：

将步骤(C12)中计算得到的大气散射系数ω₁代入天空光幅照度计算公式计算天空光幅照度E_sky1，式中：c为大气纠正前后的大气散射系数比，为测量值；

(C14)：将步骤(C12)中计算得到的大气散射系数ω₁代入大气透过率的计算公式T＝e^-ω，计算出大气透过率T₁；

(C15)：将步骤(C12)中计算得到的大气散射系数ω₁代入等式(1)右边第一项中的ω、将步骤(C14)中计算得到的大气透过率T₁代入等式(1)右边第一项中的T，计算出精度高于大气散射系数ω₁的大气散射系数ω₂；

根据步骤(C13)计算出精度更高的天空光幅照度E_sky2；

根据步骤(C14)计算出精度高于大气透过率T₁的大气透过率T₂；

依次办法，迭代计算出精度高于大气散射系数ω₂的大气散射系数ω_k、精度高于大气透过率T₂的大气透过率T_k和精度高于天空光幅照度E_sky2的天空光幅照度E_skyk；其中，k>2。

5.根据权利要求4所述的同时遥感岛礁水下地形和地物的方法，其特征在于，所述步骤(C2)中对图像中所有像元进行大气纠正的计算公式为：

$R_g=\frac{R_1-\frac{{\mathrm{\ω}}_k}{4}P\left(\mathrm{\θ}\right)}{\[T^{\′\′}+\frac{{\mathrm{\ω}}_k}{2}\]T_k}$

式中：R_g为大气纠正后得到的像元反射率，T′′＝e^-τ/cosθ,其中τ为大气光学厚度，R₁为像元表观反射率：L为像元原始辐亮度，E₀为太阳辐照度。

6.根据权利要求5所述的同时遥感岛礁水下地形和地物的方法，其特征在于，所述步骤(C3)的具体步骤为：将像元各波段的原始辐亮度减去步骤(C15)得到的天空光辐照度E_skyk与天空光水面镜面反射率R_m之积，其公式为：L₁＝L-E_skykR_m，

式中：L₁为像元消除天空光镜面反射后的辐亮度。

7.根据权利要求6所述的同时遥感岛礁水下地形和地物的方法，其特征在于，所述步骤(C4)的具体步骤为：将图像转化为反射率图像，像元所有波段的反射率均减去红外波段的反射率，其中，像元所有的波段的反射率减去红外波段的反射率的公式为：R₂＝R_g-R_ir

式中：R₂为像元消除光镜面反射后的反射率，R_ir为红外波段的反射率。

8.根据权利要求1所述的同时遥感岛礁水下地形和地物的方法，其特征在于，所述步骤(D)中水深的计算公式为：

$h=\frac{1}{\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})}\ln \[\frac{4\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})R_b-\mathrm{\β}P\left(\mathrm{\Θ}\right)}{4\mathrm{\μ}(\mathrm{\α}+\mathrm{\β})R_w-\mathrm{\β}P\left(\mathrm{\Θ}\right)}\];$

式中：P(Θ)为水体散射相函数，Θ为水体散射角；

β为水体散射系数，α为水体吸收系数，R_b为水底底质反射率，h为水深；

μ＝secθ+secφ，θ为阳光在水中传播的天顶角，φ为水中的观测天顶角；

R_w＝R_ws+R_wb，R_w为水体反射率，R_ws为水体散射光对反射率的贡献，R_wb为水底反射光对反射率的贡献。

9.根据权利要求1所述的同时遥感岛礁水下地形和地物的方法，其特征在于，步骤(D)具体为：

(D1)将步骤(B)中测量得到的某一典型地物反射率、所有波段的遥感数据代入所述物理分析模型，计算得到若干不同的水深初值；根据计算的水深初值和水深计算最佳波段规则，确定水深计算最佳波段，根据水深计算最佳波段的遥感数据算出最优水深值；

所述水深计算最佳波段规则为：当水深值在(0，0.3)米时，最佳波段为近红外波段；当水深值在[0.3，5)米时，最佳波段为近红光波段；当水深值在[5，20)米时，最佳波段为绿光波段；当水深值在[20，45]米时，最佳波段为兰光波段；

(D2)将所有波段的遥感数据、所述最优水深值代入所述物理分析模型，计算得到各波段的地物反射率，然后将计算得到的各波段的地物反射率与测量得到的地物反射率进行匹配，计算其误差；

其中，具有最小误差的地物反射率所对应的地物为该像元对应的水下地物类型，根据该地物反射率计算得到的最佳水深为该像元的水深。

10.根据权利要求9所述的同时遥感岛礁水下地形和地物的方法，其特征在于，所述步骤(D2)中计算得到的各波段的地物反射率与测量得到的地物反射率的误差的计算公式为：

${\mathrm{\ϵ}}^2=\underset{i=n}{\Σ}{(R_b-R_b^{\′})}^2;$

式中：R_b为测量得到的地物反射率，R′_b为计算得到的波段地物反射率。